Improvement of Biogas Production Utilizing a Complex Anaerobic Digestion Model and Gradient-Based Optimization

Abstract

Anaerobic digestion (AD) is a promising technology for renewable energy production from organic waste. In order to maximize the produced biogas quantity and quality, this paper deals with the optimization of the AD process in a CSTR bioreactor of a full-scale biogas plant. For this purpose, a novel approach was adopted coupling, a highly complex BioModel for AD simulation, and a gradient-based optimization method. In order to improve AD performance, the dosages of various types of biological additives, the dosages of inorganic additives, and the temperature in the bioreactor were optimized in three different scenarios. The best biogas quality was obtained using multi-objective optimization, where the objective function involves the following two conflicting objectives: the maximization of biogas production and minimization of the needed heating energy. The obtained results show that, potentially, the content of ${\mathrm{C}\mathrm{H}}_4$ can be increased by 11%, while the contents of ${\mathrm{H}}_2$, ${\mathrm{H}}_2\mathrm{S}$, and ${\mathrm{N}\mathrm{H}}_3$ can be reduced by 30%, 20%, and 81% when comparing the simulation results with the experimental data. The obtained results confirm the usefulness of the proposed approach, which can easily be adapted or upgraded for other bioreactor types.

1. Introduction

In the quest to reduce waste and produce energy, the anaerobic digestion (AD) process plays an important role [1,2]. AD is a very complex process, including hydrolysis, acidogenesis, acetogenesis, and methanogenesis, where the activity of various types of microorganisms in an oxygen-free environment degrades carbohydrates, lipids, and proteins from organic wastes into biogas [3,4]. Biogas production, as well as other important AD performances, depend on waste characteristics (waste types, composition of feeding complex substrate, chemical properties, waste pretreatment, etc.), additives (inoculum, trace elements, enzymes, etc.), and AD process conditions [5,6]. Among others, this means that by simultaneously digesting two or more organic substrates, a more favorable composition of the substrate for the AD process can be obtained. Consequently, a more stable AD process, as well as enhanced biogas and ${\mathrm{C}\mathrm{H}}_4$ production, can be achieved. Furthermore, the AD process’s parameters, such as temperature, $\mathrm{p}\mathrm{H}$ value, and additives, highly influence the formation and quality of the produced biogas and digestate [1,2]. The AD process can be performed under psychrophilic (10–25 °C), mesophilic (25–45 °C), and thermophilic (45–65 °C) temperature conditions [1,7]. Temperature affects the microorganism community, processes stability, and AD products. Some investigations show that biogas production rises with increased temperature while the content of ${\mathrm{C}\mathrm{H}}_4$ in the produced biogas reduces [3]. However, performing AD in mesophilic conditions enables the most stable AD process [3] and requires less energy to heat up the bioreactor [4]. To determine the optimal temperature for a particular AD process in a full-scale biogas plant, the interaction between temperature and other parameters, such as $\mathrm{p}\mathrm{H}$ value, nutrients, waste characteristics, and additives, should also be considered [5]. The $\mathrm{p}\mathrm{H}$ value in the bioreactor influences the growth and activities of microorganisms, the proportion of ionized and non-ionized forms of various compounds (excessive hydrogen sulfide, fatty acids, and ammonia are toxic in their non-ionized forms), and the quantity and quality of produced biogas [8]. Each enzyme and bacteria type has its own optimum range of the $\mathrm{p}\mathrm{H}$ value. For example, the optimum $\mathrm{p}\mathrm{H}$ range for acidogenic bacteria is 5.5–6.5, while the $\mathrm{p}\mathrm{H}$ from 6.5 to 8.2 is suitable for methanogenic bacteria [9]. The $\mathrm{p}\mathrm{H}$ value and its variation during the AD process depend on instantaneous ionic species concentrations in the bioreactor [10]. To assure an adequate $\mathrm{p}\mathrm{H}$ value, to overcome the inhibition of the microorganism’s activity, and to enhance the AD process, adequate additives should be used [3]. The most frequently used inorganic additives contain heavy metals, antifoaming reagents, and agents to enhance the AD process. Heavy metals may take part in various physico-chemical processes; they can form precipitates with sulfides, carbonates, and phosphates, and they can also form complexes with some AD intermediates [11]. For example, Abarghaz et al. [12] reported that a concentration at $\mathrm{C}\mathrm{a}\mathrm{C}{\mathrm{O}}_3$ of around $1.5 \mathrm{g}{\mathrm{L}}^{-1}$ is recommended for the AD process. Heavy metals, which are necessary for enzyme activation, can be stimulatory, inhibitory, or even toxic in biochemical reactions, depending on their concentrations. Many investigations show that various inorganic additives have to be added very carefully. For example, Mudhoo et al. [13] reported that a concentration of ${\mathrm{C}\mathrm{a}}^{2+}$ over $7 \mathrm{g}{\mathrm{L}}^{-1}$ can have an inhibitory effect on the AD process of synthetic waste. Guo et al. [14] showed that low concentrations of ${\mathrm{C}\mathrm{u}}^{2+}$ $\left(<0.1 \mathrm{g}{\mathrm{L}}^{-1}\right),$ ${\mathrm{F}\mathrm{e}}^{2+}\left(<4 \mathrm{g}{\mathrm{L}}^{-1}\right),$ ${\mathrm{N}\mathrm{i}}^{2+}\left(<0.1 \mathrm{g}{\mathrm{L}}^{-1}\right),$ and ${\mathrm{Z}\mathrm{n}}^{2+}\left(<0.005 \mathrm{g}{\mathrm{L}}^{-1}\right)$ promote biogas production, while higher concentrations of these heavy metals can inhibit the AD process. Furthermore, the proper addition of $\mathrm{F}\mathrm{e}{\mathrm{C}\mathrm{l}}_2$ can improve $\mathrm{C}{\mathrm{H}}_4$ content in the produced biogas [15] and reduce odor-causing ${\mathrm{H}}_2\mathrm{S}$ via the precipitation of $\mathrm{F}\mathrm{e}\mathrm{S}$ [16]. In general, various inorganic additives comprise chemical reagents, minerals, and waste sources and are able to provide micronutrients and/or promote biomass immobilization. Furthermore, the addition of biological additives, enzymes, and microorganisms, as well as the selection or manipulation of particular groups of microorganisms within the bioreactor, can improve AD performance essentially [17]. All these mentioned possibilities are investigated not only experimentally but also by numerical simulations, which provide a useful tool for AD process understanding and optimization. To better understand the AD process, reliable mechanistically inspired models for the numerical simulation of processes in an anaerobic bioreactor, based on ADM1 and BioModel, seem to be the preferred option [1,6,18,19,20,21,22,23,24]. So far, AD process optimization has been mostly performed using the genetic algorithm, particle swarm, and ant-colony algorithms [25,26,27], while gradient–based optimization methods were engaged relatively rarely [28,29]. Elagroudy et al. [28] used a gradient-based optimization method in order to find the best empirical kinetic model for the simulation of biogas production; various kinetic and performance parameters of several models (modified Gompertz, logistic model, reaction curve-type, and exponential rise to maximum) were calibrated. Furthermore, a high number of model parameters in a mechanistically inspired BioModel were calibrated by an approximation-type gradient-based optimization method quite efficiently [1,6]. A gradient-based algorithm also proved to be stable and efficient in determining AD process parameters (temperature and pH value) when using a simple mechanistically inspired model [29]. However, to the best of the authors’ knowledge, the coupling of complex mechanistically inspired mathematical models with a gradient–based algorithm in an optimization procedure of a full-scale biogas plant has not been adequately investigated and surely deserves some attention. More specifically, it is not known how a highly complex AD model would behave when engaged in gradient-based optimization. This includes unknown consequences resulting from the presumably very complicated feasible domain space of the optimization problem and the use of numerically obtained design derivatives. The novelty of this work is in trying to build such a scenario with the aim of obtaining more insight into possible difficulties and, more importantly, potential gains in developing future AD technologies.

In this paper, attention is focused on the optimization of the AD process in a full-scale biogas plant by coupling a complex mechanistically inspired BioModel and an approximation-type gradient-based optimization method. In order to improve biogas quantity and quality, three different cases of AD optimization were performed, i.e., the objective functions and constraints were defined in three ways. In all cases, the temperature and the amount of various biological and inorganic additives were used as design variables. Note that by varying inorganic additives’ dosages, one can influence the ionic species concentrations, as follows, from the corresponding acid-base equilibrium equations. In this way, the $\mathrm{p}\mathrm{H}$ value is also effectively and indirectly optimized.

2. Materials and Methods

2.1. Biogas Plant Data

In order to obtain adequate experimental data, the AD process was observed in a full-scale biogas plant in Draženci (Slovenia) for a period of one year. This plant had two mesophilic and continuously mixed reactors (CSTRs). Each single-stage CSTR had a volume of 2500 m³, and there was one common gas storage facility with a volume of 2500 m³.

The measured yearly histories of the complex substrate (F-CS) loading rate, pH value, and temperature in the CSTRs are presented in Figure 1. The $\mathrm{p}\mathrm{H}$ value throughout the AD process was $7.7\pm 0.1$, while the temperature was $42.5\pm 2.1 °\mathrm{C}$. The retention time in the CSTRs was approximately 33 days, and the AD process took place at a pressure of 1.006 bar.

The F-CS consists of the following: (i) food waste (FW) and poultry manure (PM) from a chicken farm; (ii) fat matter (FM) from a slaughterhouse; (iii) corn meal (CM) and corn silage (CS) from local farmers; and (iv) added water (W). The daily variations in FW, PM, and CS fractions in the F-CS are shown in Figure 2a for a period of one year, while the corresponding data for CM, FM, and W are presented in Figure 2b. The composition of each substrate of the F-CS,

Table 1. Composition of substrates in F-CS and the corresponding standards.

F-CS	PM	CS	CM	FM	FW	Standard SIST
Total solids, TS (%)	75.73	47.68	65.85	34.12	91.99	EN 15934:2012 [30]
Organic matter, OM (% TS)	84.76	96.58	98.35	98.00	98.08	EN 15935:2021 [31]
Carbohydrates (g/L)	10.735	6.488	6.438	4.309	31.707	ASTM E1758-01 [32]
Proteins (g/L)	45.936	33.022	56.970	74.723	66.220	ISO 937-1978 [33]
Lipids (g/L)	3.555	3.226	19.750	236.31	26.970	ISO 1443-1973 [34]
Inorganic carbon (g/L)	3.5750	2.1419	0.7153	0.1371	42.802	ASTM D4373-02 [35]
Inorganic nitrogen (g/L)	1.5532	0.6279	1.0818	1.3042	1.1523	ISO 5664:1996 [36]
Inorganic sulfur (g/L)	3.5760	0.0615	0.3745	0.0010	0.9849	ISO 9280:1990 [37]
Inorganic potassium (g/L)	3.2381	1.1058	1.1351	0.0471	2.5101	EN 16170:2016 [38]
Inorganic phosphorus (g/L)	2.6116	1.0116	3.5496	0.9006	6.6302	EN 16170:2016 [38]
Ca (g/L)	3.6139	0.1849	0.0290	1.5667	18.091	EN 16170:2016 [38]
Cr (g/L)	0.0035	0.0003	0.0008	0.0025	0.0011	EN 16170:2016 [38]
Cu (g/L)	0.0099	0.0004	0.0008	0.0024	0.0126	EN 16170:2016 [38]
Fe (g/L)	0.1103	0.0054	0.0138	0.2054	0.2752	EN 16170:2016 [38]
Mg (g/L)	0.8148	0.1084	0.3096	0.0659	0.8803	EN 16170:2016 [38]
Na (g/L)	0.3882	0.0006	0.0075	0.0618	0.5055	EN 16170:2016 [38]
NO2 (g/L)	0.0087	0.0003	0.0018	0.0001	0.0005	EN 26777:1996 [39]
Ni (g/L)	0.0005	0.0003	0.0008	0.0011	0.0011	EN 16170:2016 [38]
Pb (g/L)	0.0003	0.0003	0.0008	0.0007	0.0011	EN 16170:2016 [38]
Zn (g/L)	0.0523	0.0003	0.0054	0.0219	0.1330	EN 16170:2016 [38]

, was determined by the methods described in the corresponding standards [6]. To determine the biogas composition, the IR Gas Analyzer GA2000 plus was used.

To improve the AD process in the observed biogas plant, prevent foaming, and reduce ${\mathrm{H}}_2\mathrm{S}$ content in the produced biogas, some inorganic (SensoPower Liquid, purchased from Phytobiotics Futterzusatzstoffe GmbH, Eltville, Germany and Kemira BDP-840, purchased from Kemira Oyj, Helsinki, Finland) and biological (SensoPower Flex, purchased from Phytobiotics Futterzusatzstoffe GmbH, Eltville, Germany) additives were daily added to the F-CS. The SensoPower Liquid (SPL) contains all essential trace elements ($\mathrm{N}\mathrm{i}$, $\mathrm{C}\mathrm{o}$, $\mathrm{N}\mathrm{a}$, $\mathrm{C}\mathrm{u}$, and $\mathrm{S}\mathrm{e}$), which are important for the healthy development of the microorganisms; furthermore, they enable efficient fermentation and prevent foam development. The Kemira BDP-840 (Kemira) contains $\mathrm{F}\mathrm{e}{\mathrm{C}\mathrm{l}}_2$ which reduces the production of ${\mathrm{H}}_2\mathrm{S}$ during the AD process. The SensoPower Flex (SPF) contains enzymes (cellulase, xylanase, endo-1,4, …) to enhance hydrolysis in the AD process. The amounts of the added SPL, Kemira, and SPF are $1.0 \mathrm{L} {\mathrm{d}\mathrm{a}\mathrm{y}}^{-1}$, $100 \mathrm{L} {\mathrm{d}\mathrm{a}\mathrm{y}}^{-1}$, and $1500 \mathrm{g} {\mathrm{d}\mathrm{a}\mathrm{y}}^{-1}$, respectively.

The measured total unrefined biogas volume in one year was about 4,370,000 m³. The produced biogas contained approximately 54% ${\mathrm{C}\mathrm{H}}_4$, 45% ${\mathrm{C}\mathrm{O}}_2$, 60 ppm of ${\mathrm{H}}_2$, 200 ppm of ${\mathrm{H}}_2\mathrm{S}$, and 800 ppm of ${\mathrm{N}\mathrm{H}}_3$.

2.2. Numerical Simulation of the AD Process

To simulate the AD process, an in-house developed complex BioModel, which is described in detail in [1,6], was used. To summarize briefly, this model takes into account chemical (acid-base reactions), biochemical (hydrolysis, acidogenesis, acetogenesis, and methanogenis), and physico-chemical (liquid-gas and liquid-solid mass transfers) processes, where carbohydrates, proteins, and lipids degrade into ${\mathrm{C}\mathrm{H}}_4$ and many other by-products. The model contains 184 model parameters, which were calibrated by an active set optimization procedure with respect to the experimental data obtained from a full-scale biogas plant [1,6]. In this model, the $\mathrm{p}\mathrm{H}$ value was treated as a state variable and determined by the corresponding algebraic and differential equations.

The general assumptions are as follows: (i) biogas, containing $\mathrm{C}{\mathrm{H}}_4, \mathrm{C}{\mathrm{O}}_2, {\mathrm{H}}_2, {\mathrm{H}}_2\mathrm{S},\mathrm{N}{\mathrm{H}}_3$, and ${\mathrm{H}}_2\mathrm{O}$ vapor, follows the ideal gas behavior; (ii) system pressure is constant; (iii) ${\mathrm{H}}_2\mathrm{O}$ vapor in biogas is at a saturation state; (iv) extracellular biochemical reactions are based on modified Michaelis–Menten kinetics; (v) intracellular biochemical reactions are based on Monod kinetics, modified including non-competitive microbial growth inhibitions; (vi) acid-base recombinations occur instantaneously; (vii) the amounts of added extracellular enzymes, which affect the hydrolysis of carbohydrates, proteins, and lipids, are equal; and (viii) the characteristics of a particular feedstock substrate are constant.

2.3. Optimization of the AD Process

In order to achieve the best AD performances of the considered biogas plant with prescribed feeding substrate composition and feeding strategy during a time period of 365 days, the following AD process parameters can be optimized: (i) biological additives (initial concentrations of bacteria groups and dosages of enzymes), (ii) inorganic additives, which are involved in the physico-chemical processes and affect the $\mathrm{p}\mathrm{H}$ value; and (iii) temperature in the bioreactor.

The considered optimization problem can be defined as follows:

$$\min g_0(\mathbf{x},\mathbf{q})$$

(1)

$$g_i(\mathbf{x},\mathbf{q})\le 0, i=1,\dots ,k^{*}$$

(2)

$$h\left(\mathbf{x},\mathbf{q},\dot{\mathbf{q},}t\right)=0$$

(3)

where the vectors $\mathbf{x}\mathsf{ϵ}{ R}^{n^{*}}$ and $\mathbf{q}\mathsf{ϵ}{ R}^{m^{*}}$ assemble the design and state variables, respectively; the time derivatives of the state variables are denoted by $\dot{\mathbf{q}}\mathsf{ϵ}{ R}^{m^{*}}$. The dependency of $\mathbf{q}$ on $t$ and $\mathbf{x}$ is given by the state equation, Equation (3), which is defined by the whole BioModel. The symbols $g_0$ and $g_i$ denote the objective and constraint functions, respectively. These functions involve various AD performance quantities and process parameters. The symbols $n^{*}$, $k^{*}$, and $m^{*}$ denote the number of design variables, constraints, and state variables, respectively.

The design variables $x_i$ $i={1\dots n}^{*}$ ($n^{*}=$17) are related to the following:

• Biological additives:

  √

  The initial concentrations of bacteria groups are $X_{i,\mathrm{i}\mathrm{n}}, i\in I_{\mathrm{b}\mathrm{a}\mathrm{c}}, I_{\mathrm{b}\mathrm{a}\mathrm{c}}=$ {Asu, Aaa, Agly, Aoa, Apro, Abu, Ava, Mac, Mhyd, Ss, Spro, Sac, Shyd}; these bacteria groups are acidogenic degraders of sugar (Asu), amino acids (Aaa), glycerol (Agly), and oleic acid (Aoa), acetogenic degraders of propionic acid (Apro), butyric acid (Abu), and valeric acid (Ava), methanogenic degraders bacteria of acetate (Mac), hydrogen (Mhyd), and sulfate-reducing bacteria involved in the reduction in sulfate (Ss), propionate (Spro), acetate (Sac), and hydrogen (Shyd),

  √

  The dosage of the SensoPower Flex additive $c_{\mathrm{A}\mathrm{d}\mathrm{d}\mathrm{S}\mathrm{P}\mathrm{F}}$, which contains enzymes;

• Inorganic additives;

  √

  The dosage of the Kemira BDP-840 additive $c_{\mathrm{A}\mathrm{d}\mathrm{d}\mathrm{K}\mathrm{e}\mathrm{m}},$ which contains $\mathrm{F}\mathrm{e}{\mathrm{C}\mathrm{l}}_2$;

  √

  The dosage of the SensoPower Liquid additive $c_{\mathrm{A}\mathrm{d}\mathrm{d}\mathrm{S}\mathrm{P}\mathrm{L}}$, which contains trace elements ($\mathrm{N}\mathrm{i}$, $\mathrm{C}\mathrm{o}$, $\mathrm{N}\mathrm{a}$, $\mathrm{C}\mathrm{u}$, and $\mathrm{S}\mathrm{e}$);

• Temperature $T$.

To gain a better insight into the AD optimization process, three different optimization tasks were defined; these tasks are denoted in the following as cases A, B, and C.

Case A. The goal in Case A was to maximize the biogas production. Since the initial concentration of the complete microbial consortium has to be limited, the maximal concentration of all microbial groups in the feeding substrate was constrained. Therefore, the objective function $g_0$ and one constraint function $g_1$ are defined as follows:

$$g_0=-\frac{1}{t_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}-t_{\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{b}}}\frac{{\int }_{t_{\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{b}}}^{t_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}}{Q}_{\mathrm{b}\mathrm{i}\mathrm{o}\mathrm{g}\mathrm{a}\mathrm{s}}\left(t\right)dt}{ K_{\mathrm{b}\mathrm{i}\mathrm{o}\mathrm{g}\mathrm{a}\mathrm{s}}}$$

(4)

$$g_1=\frac{\sum _i{X}_i-X_{\mathrm{b}\mathrm{a}\mathrm{c}}^{\mathrm{m}\mathrm{a}\mathrm{x}}}{X_{\mathrm{b}\mathrm{a}\mathrm{c}}^{\mathrm{m}\mathrm{a}\mathrm{x}}}, i\in I_{\mathrm{b}\mathrm{a}\mathrm{c}}$$

(5)

where $Q_{\mathrm{b}\mathrm{i}\mathrm{o}\mathrm{g}\mathrm{a}\mathrm{s}}$ denotes the biogas flow rate, $K_{\mathrm{b}\mathrm{i}\mathrm{o}\mathrm{g}\mathrm{a}\mathrm{s}}$ is the corresponding normalization constant, $t_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}$ represents the total time of the AD process observation, $t_{\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{b}}$ is the AD process stabilization time, $X_i$ is the initial concentration of the $i$th bacteria group, $i\in I_{\mathrm{b}\mathrm{a}\mathrm{c}}$, while $X_{\mathrm{b}\mathrm{a}\mathrm{c}}^{\mathrm{m}\mathrm{a}x}$ represents the maximal initial concentration of the sum of all bacteria groups.

Case B. Since the presence of ${\mathrm{H}}_2\mathrm{S}$ in the biogas produced by AD can be high (up to 30,000 ppm), it can lead to many problems in biogas storage and its further use (boilers, turbines, gas engines, dual-fuel diesel engines, full cells…) [4]. Note that besides causing the corrosion of various components, ${\mathrm{H}}_2\mathrm{S}$ is also toxic to methanogenic microorganisms. Therefore, the content of ${\mathrm{H}}_2\mathrm{S}$ in the produced biogas has to be limited. The goal in Case B is to maximize the biogas production while imposing the limitation of the total concentration of all bacteria groups and the content of hydrogen sulfide in the produced biogas; the corresponding objective and constraint functions $g_0$ and $g_1$ are given by Equations (4) and (5), while the constraint function $g_2$ is defined as follows:

$$g_2=\frac{{\int }_{t_{\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{b}}}^{t_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}}{Q}_{{\mathrm{H}}_2\mathrm{S}}\left(t\right)dt}{{\int }_{t_{\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{b}}}^{t_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}}{Q}_{\mathrm{b}\mathrm{i}\mathrm{o}\mathrm{g}\mathrm{a}\mathrm{s}}\left(t\right)dt}-{\phi }_{{\mathrm{H}}_2\mathrm{S},\mathrm{m}\mathrm{a}\mathrm{x}}$$

(6)

where $Q_{{\mathrm{H}}_2\mathrm{S}}$ denotes the flow rate of ${\mathrm{H}}_2\mathrm{S}$, while ${\phi }_{{\mathrm{H}}_2\mathrm{S},\mathrm{m}\mathrm{a}\mathrm{x}}$ is the prescribed maximal content of ${\mathrm{H}}_2\mathrm{S}$ in the produced biogas.

Case C. Achieving an adequate AD process temperature for maximal biogas production would typically require heating up the bioreactors substantially, which may not be economically justified [4]. Furthermore, it is worth noting that low AD process temperatures improve the process stability and inactivation of the pathogens. Since the cost of heating the bioreactor is desired to be minimal, this needs to be reflected in the definition of the optimization problem. This means that there are the following two objectives: biogas production maximization and heating cost minimization. Consequently, a multi-objective optimization problem has to be addressed. This is usually conducted by solving an adequate single-objective problem, which is obtained by summing individual weighted objective functions. In this work, the single-objective function was defined as follows:

$$g_0={-\psi }_1\frac{1}{t_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}-t_{\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{b}}}\underset{t_{\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{b}}}{\overset{t_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}}{\int }}\frac{Q_{\mathrm{b}\mathrm{i}\mathrm{o}\mathrm{g}\mathrm{a}\mathrm{s}}\left(t\right)}{K_{\mathrm{b}\mathrm{i}\mathrm{o}\mathrm{g}\mathrm{a}\mathrm{s}}}dt+{\psi }_2\frac{T}{K_{\mathrm{T}}}$$

(7)

where ${\psi }_1$ and ${\psi }_2$ are the weighing factors for biogas production and heating cost, $T$ is the temperature of the AD process, and $K_{\mathrm{T}}$ is the temperature normalization constant. The objective in Case C was to maximize the biogas flow rate and minimize the bioreactor heating energy, while satisfying seven constraints. The first constraint was related to the total concentration of all bacteria groups; Equation (5). The next three constraint functions, which were related to the maximum allowed fractions of ${\mathrm{H}}_2\mathrm{S}$, ${\mathrm{H}}_2$, and ${\mathrm{N}\mathrm{H}}_3$ in the produced biogas, were defined by Equation (6) for ${\mathrm{H}}_2\mathrm{S}$, Equation (8) for ${\mathrm{H}}_2$, and Equation (9) for ${\mathrm{N}\mathrm{H}}_3$. Due to the fact that the amount of ${\mathrm{C}\mathrm{H}}_4$ should represent the highest part of the produced biogas, the minimal allowed fraction of ${\mathrm{C}\mathrm{H}}_4$ in the produced biogas was involved into the constraint function, which is written as Equation (10). Finally, the $\mathrm{p}\mathrm{H}$ values, which are dynamically computed in dependence of various physico-chemical processes during the AD process [1,6], are limited to be within the allowed interval by Equations (11) and (12). The constraints (3)–(7) were defined as follows:

$$g_3=\frac{{\int }_{t_{\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{b}}}^{t_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}}{Q}_{{\mathrm{H}}_2}\left(t\right)dt}{{\int }_{t_{\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{b}}}^{t_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}}{Q}_{\mathrm{b}\mathrm{i}\mathrm{o}\mathrm{g}\mathrm{a}\mathrm{s}}\left(t\right)dt}-{\phi }_{{\mathrm{H}}_2,\mathrm{m}\mathrm{a}\mathrm{x}}$$

(8)

$$g_4=\frac{{\int }_{t_{\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{b}}}^{t_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}}{Q}_{{\mathrm{N}\mathrm{H}}_3}\left(t\right)dt}{{\int }_{t_{\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{b}}}^{t_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}}{Q}_{\mathrm{b}\mathrm{i}\mathrm{o}\mathrm{g}\mathrm{a}\mathrm{s}}\left(t\right)dt}-{\phi }_{{\mathrm{N}\mathrm{H}}_3,\mathrm{m}\mathrm{a}\mathrm{x}}$$

(9)

$$g_5={\phi }_{{\mathrm{C}\mathrm{H}}_4,\mathrm{m}\mathrm{i}\mathrm{n}}-\frac{{\int }_{t_{\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{b}}}^{t_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}}{Q}_{{\mathrm{C}\mathrm{H}}_4}\left(t\right)dt}{{\int }_{t_{\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{b}}}^{t_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}}{Q}_{\mathrm{b}\mathrm{i}\mathrm{o}\mathrm{g}\mathrm{a}\mathrm{s}}\left(t\right)dt}$$

(10)

$$g_6=\frac{{\mathrm{p}\mathrm{H}}_{\mathrm{m}\mathrm{i}\mathrm{n}}-\frac{{\int }_{t_{\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{b}}}^{t_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}}\mathrm{p}\mathrm{H}\left(t\right)dt}{t_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}-t_{\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{b}}}}{{\mathrm{p}\mathrm{H}}_{\mathrm{m}\mathrm{i}\mathrm{n}}}$$

(11)

$$g_7=\frac{\frac{{\int }_{t_{\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{b}}}^{t_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}}\mathrm{p}\mathrm{H}\left(t\right)dt}{t_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}-t_{\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{b}}}-{\mathrm{p}\mathrm{H}}_{\mathrm{m}\mathrm{a}\mathrm{x}}}{ {\mathrm{p}\mathrm{H}}_{\mathrm{m}\mathrm{a}\mathrm{x}}}$$

(12)

where $Q_{{\mathrm{H}}_2}$, $Q_{{\mathrm{N}\mathrm{H}}_3}$, and $Q_{{\mathrm{C}\mathrm{H}}_4}$ denote the flow rates of ${\mathrm{H}}_2$, ${\mathrm{N}\mathrm{H}}_3$, and ${\mathrm{C}\mathrm{H}}_4$, respectively; the symbols ${\phi }_{{\mathrm{H}}_2,\mathrm{m}\mathrm{a}\mathrm{x}}$ and ${\phi }_{{\mathrm{N}\mathrm{H}}_3,\mathrm{m}\mathrm{a}\mathrm{x}}$ denote the prescribed maximal contents of ${\mathrm{H}}_2\mathrm{S}$ and ${\mathrm{N}\mathrm{H}}_3$ in biogas, while ${\phi }_{{\mathrm{C}\mathrm{H}}_4,\mathrm{m}\mathrm{i}\mathrm{n}}$ represents the minimal allowed part of ${\mathrm{C}\mathrm{H}}_4$ in the produced biogas; the symbols ${\mathrm{p}\mathrm{H}}_{\mathrm{m}\mathrm{i}\mathrm{n}}$ and ${\mathrm{p}\mathrm{H}}_{\mathrm{m}\mathrm{a}\mathrm{x}}$ denote the prescribed minimal and maximal $\mathrm{p}\mathrm{H}$ values.

3. Results and Discussion

The BioModel and the optimization procedure were developed in-house and coded in the C# computer language. The system of differential equations was integrated by the Runge–Kutta method. The utilized optimizer was based on an approximation method [40], which generated a sequence of approximate optimization problems. These approximate problems are strictly convex and separable and can be solved rather efficiently for a sequence of approximate solutions. Since the conservativeness of the approximate objective and constraint functions is adjusted automatically, the solution to the process is usually fast and stable. A gradient-based optimizer clearly needs the design derivatives of the involved functions. Since analytical derivatives cannot be obtained easily and any further AD model updates would invalidate the derived formulas, numerical differentiation using forward differences was used in this work. Numerical differentiation can be computationally quite costly. Therefore, the computation of design derivatives was parallelized. By utilizing parallelization, one full optimization cycle took about 1 min of the processing time on an i7 CPU computer.

3.1. Numerical Simulation of the AD Process

The numerical simulation of the AD process was performed using the BioModel, which was thoroughly calibrated and validated in our previous work [1,6]. Just to illustrate the performance of the model, a few comparisons of simulated and experimental data are shown in the following. Figure 3 shows the comparison of biogas and ${\mathrm{C}\mathrm{H}}_4$ flow rates, while the pH value as well as ${\mathrm{H}}_2\mathrm{S}$ and ${\mathrm{H}}_2$ flow rates are shown in Figure 4. The differences between simulated and experimental average values (averaged over the whole period of the AD process) of biogas, ${\mathrm{C}\mathrm{H}}_4$, ${\mathrm{H}}_2$, and ${\mathrm{H}}_2\mathrm{S}$ flow rates are all in the range between 0.6% and 2.2%. Furthermore, the average p$\mathrm{H}$ value differs negligibly, that is, less than 0.2%.

Note that the AD model used in this work was also rigorously validated [1,6]. For example, during validation, the simulated biogas production rate per unit of the bioreactor volume $\left(2.31\pm 0.25 \mathrm{L}{\mathrm{L}}^{-1} {\mathrm{d}\mathrm{a}\mathrm{y}}^{-1}\right),$ $\mathrm{p}\mathrm{H}$ value $\left(7.70\pm 0.02\right)$, the ${\mathrm{C}\mathrm{H}}_4$ content in biogas $\left(54.4\pm 0.88\mathrm{\%}\right)$, ${\mathrm{H}}_2$ content $\left(62.11\pm 5.55 \mathrm{p}\mathrm{p}\mathrm{m}\right)$, and ${\mathrm{H}}_2\mathrm{S}$ content $\left(195.89\pm 19.92 \mathrm{p}\mathrm{p}\mathrm{m}\right)$ were very close to the measured values. All statistical indicators (the mean absolute error ${\epsilon }_{\mathrm{M}\mathrm{A}\mathrm{E}}$, root mean square error ${\epsilon }_{\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}$, coefficient of determination $R^2$, and the relative index of agreement $I_{\mathrm{A},\mathrm{r}\mathrm{e}\mathrm{l}}$) were satisfactory; for example, the relative index of agreement $I_{\mathrm{A},\mathrm{r}\mathrm{e}\mathrm{l}}$ was about 0.95 for biogas and ${\mathrm{C}\mathrm{H}}_4$ flow rates and about 0.75 for ${\mathrm{H}}_2$ and ${\mathrm{H}}_2\mathrm{S}$ flow rates.

According to the above discussion, one can reasonably assume that the highly complex BioModel, which involves 184 well-calibrated model parameters, represents a promising foundation for AD process optimization.

3.2. Optimization of the AD Process

The numerical optimization process was carefully observed in all considered cases. Numerical instabilities were not observed, and the optimization history was relatively smooth. Furthermore, there were no indications pointing to possible problems related to using finite differences in computing the design derivatives. All considered cases were solved successfully; the number of optimization cycles needed to obtain the optimum values of design variables was 63 in Case A, 47 in Case B, and 45 in Case C.

The prescribed values of limit parameters in the constraint functions were as follows: $X_{\mathrm{b}\mathrm{a}\mathrm{c}}^{\mathrm{m}\mathrm{a}x}=5 \mathrm{g}{\mathrm{L}}^{-1}$ (Cases A, B, and C), ${\phi }_{{\mathrm{H}}_2\mathrm{S},\mathrm{m}\mathrm{a}\mathrm{x}}=200 \mathrm{p}\mathrm{p}\mathrm{m}$ (Case B) and ${\phi }_{{\mathrm{H}}_2\mathrm{S},\mathrm{m}\mathrm{a}\mathrm{x}}=150 \mathrm{p}\mathrm{p}\mathrm{m}$ (Case C), ${\phi }_{{\mathrm{H}}_2,\mathrm{m}\mathrm{a}\mathrm{x}}=50 \mathrm{p}\mathrm{p}\mathrm{m}$ (Case C), ${\phi }_{{\mathrm{N}\mathrm{H}}_3,\mathrm{m}\mathrm{a}\mathrm{x}}=500 \mathrm{p}\mathrm{p}\mathrm{m}$ (Case C), ${\phi }_{{\mathrm{C}\mathrm{H}}_4,\mathrm{m}\mathrm{i}\mathrm{n}}=60\%,$ ${\mathrm{p}\mathrm{H}}_{\mathrm{m}\mathrm{i}\mathrm{n}}=5.5$ (Case C)$,$ and ${\mathrm{p}\mathrm{H}}_{\mathrm{m}\mathrm{a}\mathrm{x}}=8.5$ (Case C). The normalization constants were $K_{\mathrm{b}\mathrm{i}\mathrm{o}\mathrm{g}\mathrm{a}\mathrm{s}}={10}^{-8}$ and $K_{\mathrm{T}}=100$ while the total simulation and stabilization times were $t_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}=500 \mathrm{d}\mathrm{a}\mathrm{y}\mathrm{s}$ and $t_{\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{b}}=135 \mathrm{d}\mathrm{a}\mathrm{y}\mathrm{s}$.

During optimization, the values of design variables $x_i$ were allowed to vary between the prescribed lower $x_i^{\mathrm{l}\mathrm{o}}$ and upper $x_i^{\mathrm{u}\mathrm{p}}$ limits. These limits, initial values, and optimal values for Cases A, B, and C are collected in

Table 2. Lower and upper limits, initial, and optimal values of design variables.

No	$\mathbf{Design} \mathbf{Variable} {\mathit{x}}_{\mathit{i}}$	${\mathit{x}}_{\mathit{i}}^{\mathbf{l}\mathbf{o}}$	${\mathit{x}}_{\mathit{i}}^{\mathbf{u}\mathbf{p}}$	Initial State	Optimal State
					Case A	Case B	Case C
1	$X_{\mathrm{A}\mathrm{s}\mathrm{u},\mathrm{i}\mathrm{n}}\left(\mathrm{g}{\mathrm{L}}^{-1}\right)$	0.1	0.5	0.2806	0.3077	0.2510	0.2876
2	$X_{\mathrm{A}\mathrm{a}\mathrm{a},\mathrm{i}\mathrm{n}}\left(\mathrm{g}{\mathrm{L}}^{-1}\right)$	0.1	0.5	0.2802	0.3057	0.2493	0.2864
3	$X_{\mathrm{A}\mathrm{g}\mathrm{l}\mathrm{y},\mathrm{i}\mathrm{n}}\left(\mathrm{g}{\mathrm{L}}^{-1}\right)$	0.1	0.5	0.2804	0.3073	0.2508	0.2874
4	$X_{\mathrm{A}\mathrm{o}\mathrm{a},\mathrm{i}\mathrm{n}}\left(\mathrm{g}{\mathrm{L}}^{-1}\right)$	0.1	0.5	0.2808	0.3108	0.2536	0.2890
5	$X_{\mathrm{A}\mathrm{p}\mathrm{r}\mathrm{o},\mathrm{i}\mathrm{n}}\left(\mathrm{g}{\mathrm{L}}^{-1}\right)$	0.1	0.5	0.2891	0.3410	0.2966	0.3249
6	$X_{\mathrm{A}\mathrm{b}\mathrm{u},\mathrm{i}\mathrm{n}}\left(\mathrm{g}{\mathrm{L}}^{-1}\right)$	0.1	0.5	0.2916	0.3552	0.3158	0.3382
7	$X_{\mathrm{A}\mathrm{v}\mathrm{a},\mathrm{i}\mathrm{n}}\left(\mathrm{g}{\mathrm{L}}^{-1}\right)$	0.1	0.5	0.2881	0.3716	0.3355	0.3498
8	$X_{\mathrm{M}\mathrm{a}\mathrm{c},\mathrm{i}\mathrm{n}}\left(\mathrm{g}{\mathrm{L}}^{-1}\right)$	0.1	0.5	0.2652	0.4948	0.4680	0.3385
9	$X_{\mathrm{M}\mathrm{h}\mathrm{y}\mathrm{d},\mathrm{i}\mathrm{n}}\left(\mathrm{g}{\mathrm{L}}^{-1}\right)$	0.1	0.5	0.1012	0.1264	0.1000	0.1096
10	$X_{\mathrm{S}\mathrm{p}\mathrm{r}\mathrm{o},\mathrm{i}\mathrm{n}}\left(\mathrm{g}{\mathrm{L}}^{-1}\right)$	0.1	0.5	0.2438	0.2546	0.1915	0.2468
11	$X_{\mathrm{S}\mathrm{a}\mathrm{c},\mathrm{i}\mathrm{n}}\left(\mathrm{g}{\mathrm{L}}^{-1}\right)$	0.1	0.5	0.2801	0.3021	0.2482	0.2863
12	$X_{\mathrm{S}\mathrm{s},\mathrm{i}\mathrm{n}}\left(\mathrm{g}{\mathrm{L}}^{-1}\right)$	0.1	0.5	0.2805	0.3036	0.2493	0.2867
13	$X_{\mathrm{S}\mathrm{h}\mathrm{y}\mathrm{d},\mathrm{i}\mathrm{n}}\left(\mathrm{g}{\mathrm{L}}^{-1}\right)$	0.1	0.5	0.2671	0.4712	0.2347	0.2730
14	$c_{\mathrm{A}\mathrm{d}\mathrm{d}\mathrm{S}\mathrm{P}\mathrm{F}}\left(\mathrm{g}{\mathrm{L}}^{-1}\right)$	0.001	0.1	0.0060	0.0678	0.0517	0.0453
15	$c_{\mathrm{A}\mathrm{d}\mathrm{d}\mathrm{K}\mathrm{e}\mathrm{m}}\left(\mathrm{g}{\mathrm{L}}^{-1}\right)$	0.1	5.0	0.4500	4.9999	4.9863	4.9307
16	$c_{\mathrm{A}\mathrm{d}\mathrm{d}\mathrm{S}\mathrm{P}\mathrm{L}}\left(\mathrm{g}{\mathrm{L}}^{-1}\right)$	0.001	0.1	0.0050	0.0140	0.0010	0.0076
17	$T(°\mathrm{C})$	25.0	65.0	42.100	59.024	45.288	25.001

. The initial values of the concentrations of all microbial groups were determined previously by the calibration of BioModel parameters [1,6]. The initial values of the amounts of additives were taken from the actual data of the considered biogas plant. The initial value of temperature was taken as the average temperature in the bioreactor in an observed time period of 365 days. The AD performances obtained with the initial values of design variables (initial state) agree very well with the measured AD performances presented in Figure 3 and Figure 4.

With respect to the initial state, the obtained optimal values of the initial concentration of all bacteria were higher by 18%, 0.5%, and 8% in Cases A, B, and C, respectively. Furthermore, optimization delivered the highest differences in the initial concentration of methanogenic acetate degraders ($X_{\mathrm{M}\mathrm{a}\mathrm{c},\mathrm{i}\mathrm{n}}$), followed by the acetogenic valerate degraders ($X_{\mathrm{A}\mathrm{v}\mathrm{a},\mathrm{i}\mathrm{n}}$). With respect to the initial state, the optimized $X_{\mathrm{M}\mathrm{a}\mathrm{c},\mathrm{i}\mathrm{n}}$ were higher by 87%, 77%, and 28% in Cases A, B, and C, respectively, while the optimized $X_{\mathrm{A}\mathrm{v}\mathrm{a},\mathrm{i}\mathrm{n}}$ in Cases A, B, and C were higher by 29%, 16%, and 21%, respectively. The highest concentrations were obtained for the methanogenic acetate degraders in Cases A and B, where the objective function is related to biogas production only.

The obtained optimal bacteria composition of biological additive is presented in Figure 5, where the relative abundance of each bacterium as well as each bacteria group is presented. By comparing the relative abundance of each bacterium in the optimal states in Cases A, B, and C with the initial state, the highest differences were observed for methanogenic acetate degraders. The parts of individual bacteria in the acidogens, acetogens, and SRB groups varied within about 5% in all optimization cases. The dominant bacteria in the methanogenic group are methanogenic degraders of acetate ($\mathrm{M}\mathrm{a}\mathrm{c}$), which at an initial state represent 72% of methanogenic bacteria. After optimization, the methanogenic degraders of acetate represent 80%, 82%, and 76% of methanogenic bacteria in Cases A, B, and C, respectively. By comparing the relative abundance of each bacteria group after the optimization, it was evident that the parts of all acidogens and SRB bacteria in all cases were lower than at their initial state, while the parts of all methanogens and acetogens were higher with respect to the initial state. An increase in methanogens and acetogens was expected since methanogens play a key role in biogas production during the AD process, while acetogens are essential to ensure process stability [41].

The optimal dosages of the Kemira BDP-840 additive were, in all cases, higher than at the initial state (experimental data). This additive promotes biogas production (Case A) as well as lower ${\mathrm{H}}_2\mathrm{S}$ content in the produced biogas (Cases B and C). The optimized dosage of the enzymes provided by the SPF additive is, in all cases, higher than that actually added in the real biogas plant, Figure 6. The optimal dosage of the inorganic additive SPL was the lowest in Case B, followed by Cases C and A. Because of the lowest optimized values of the amount of SPL, it can be concluded that the conditions in the bioreactor in Case B are the most beneficial for the AD process; note that higher concentrations of inorganic additives can inhibit the AD process [14].

To estimate potential AD performance improvements, the quantity and quality of the produced biogas were analyzed in all cases by considering the impact of temperature, biological and inorganic additives, and their interactions. Some results with the most important findings are presented in the following.

In Cases A, B, and C, AD performances obtained by the optimized values of design variables were compared with those obtained by their initial values (initial state). The biogas flow rate for initial and optimal states is presented in Figure 7. The highest biogas flow rate reached Case A, where it increased by 30%. Meanwhile, in Cases B and C, the biogas production increased by about 14% and 1%, respectively. In Case A, the optimized temperature was increased to 60 $°\mathrm{C}$. This high temperature is the main reason for increased biogas production, which is in accordance with the fact that thermophilic temperature can promote biogas production [4]. The high optimized dosage $c_{\mathrm{A}\mathrm{d}\mathrm{d}\mathrm{K}\mathrm{e}\mathrm{m}}$ of the Kemira BDP-840 additive is another important reason for the highest biogas flow rate in Case A [42]. In Case B, where ${\mathrm{H}}_2\mathrm{S}$ production was limited, the obtained optimal value of temperature was around 45 $°\mathrm{C}$. In Case C, where the objective was to increase biogas production and reduce the bioreactor heating energy, the obtained optimal temperature was about 25 $°\mathrm{C}$; consequently, the biogas production rate could not increase essentially with respect to the initial state. Anyhow, by optimizing the dosage of biological and inorganic additives, it was possible to preserve the initial biogas production rate while ensuring a substantial decrease in temperature (lower heating cost). The relative biogas production rate (the production rate with respect to bioreactor volume) was $2.41 \mathrm{L}{\mathrm{L}}^{-1} {\mathrm{d}\mathrm{a}\mathrm{y}}^{-1}$ in Case C, while the highest value of $3.10 \mathrm{L}{\mathrm{L}}^{-1} {\mathrm{d}\mathrm{a}\mathrm{y}}^{-1}$ was obtained in Case A.

The ${\mathrm{C}\mathrm{H}}_4$ flow rate and ${\mathrm{C}\mathrm{H}}_4$ content in biogas for initial and optimal states are presented in Figure 8. The highest ${\mathrm{C}\mathrm{H}}_4$ production rate and content in biogas were reached in Case C. With respect to the initial state, the ${\mathrm{C}\mathrm{H}}_4$ production rate was higher by 9%, 11%, and 12.5% in Cases A, B, and C, respectively. Figure 8 shows that there was no significant difference in ${\mathrm{C}\mathrm{H}}_4$ flow rates among all three cases, but the content of ${\mathrm{C}\mathrm{H}}_4$ in the produced biogas increased from Case A through Case B to Case C. The obtained results are in agreement with those presented by Cao et al. [43], who demonstrated that the highest content of ${\mathrm{C}\mathrm{H}}_4$ could be reached at mesophilic temperature conditions. The best relative ${\mathrm{C}\mathrm{H}}_4$ production rate of $1.4 \mathrm{L}{\mathrm{L}}^{-1} {\mathrm{d}\mathrm{a}\mathrm{y}}^{-1}$ was reached in Case C.

The ${\mathrm{H}}_2$ flow rate and ${\mathrm{H}}_2$ content in biogas for initial and optimal states are presented in Figure 9. The smallest ${\mathrm{H}}_2$ production was reached in Case C, where the relative ${\mathrm{H}}_2$ production rate was lower by about 29% with respect to the initial state. With respect to the initial state, the relative ${\mathrm{H}}_2$ production rate increased by about 11% in Case B and by about 87% in Case A. The highest relative ${\mathrm{H}}_2$ production rate and the highest ${\mathrm{H}}_2$ content in the produced biogas are obtained in Case A, where the optimized value of the temperature was the highest, and the optimized additives concentrations led to a low pH value. The obtained results are in agreement with Wang et al. [44], who reported that a thermophilic condition and low pH value lead to a higher hydrogen production rate.

The ${\mathrm{H}}_2\mathrm{S}$ flow rate and ${\mathrm{H}}_2\mathrm{S}$ content in biogas at initial and optimal states are presented in Figure 10. The smallest ${\mathrm{H}}_2\mathrm{S}$ production was reached in Case C. With respect to the initial state, the relative ${\mathrm{H}}_2\mathrm{S}$ production rate was higher by about 114% in Case A and by about 22% in Case B. Only in Case C was the obtained relative ${\mathrm{H}}_2\mathrm{S}$ production rate lower by about 20% when compared to the initial state. The lowest ${\mathrm{H}}_2\mathrm{S}$ content in the produced biogas in Case C is the consequence of the interactions of low temperature and appropriate amounts of biological and inorganic additives, which influenced the $\mathrm{p}\mathrm{H}$ value favorably.

The ${\mathrm{N}\mathrm{H}}_3$ flow rate and the ${\mathrm{N}\mathrm{H}}_3$ content in biogas for initial and optimal states are presented in Figure 11. The smallest ${\mathrm{N}\mathrm{H}}_3$ production was reached in Case C. The relative ${\mathrm{N}\mathrm{H}}_3$ production rate increased by about 270% in Case A and by about 43% in Case B with respect to the initial state; only Case C delivered the relative production rate of ${\mathrm{N}\mathrm{H}}_3,$ which was about 80% lower than at the initial state.

A comparison of the simulated $\mathrm{p}\mathrm{H}$ value during the AD process at initial and optimal states in Cases A, B, C is presented in Figure 12. The average pH value during the 365 days of the AD process was 7.7 for the initial design; this value is close to the experimental one. The optimized average $\mathrm{p}\mathrm{H}$ values were 7.5 in Case A, 7.6 in Case B, and 7.7 in Case C.

The variations in the production of biogas and its compounds at optimal states (Cases A, B, and C) with respect to the initial state are presented in Figure 13. It can be seen that the increased biogas production in Case A is a consequence of the increase in all biogas compounds: ${\mathrm{C}\mathrm{H}}_4$, ${\mathrm{H}}_2$, ${\mathrm{H}}_2\mathrm{S}$, ${\mathrm{N}\mathrm{H}}_3$, ${\mathrm{C}\mathrm{O}}_2$, and ${\mathrm{H}}_2\mathrm{O}$ vapour. The main part of the produced biogas is ${\mathrm{C}\mathrm{O}}_2$. Furthermore, in Case A, the content of ${\mathrm{H}}_2\mathrm{O}$ vapour is the highest (more than 10%), while in Case C, it represents only 3%. This is expected since biogas typically contains 10% by volume of ${\mathrm{H}}_2\mathrm{O}$ vapour at 43 $°\mathrm{C}$, 5% by volume at 32 $°\mathrm{C}$, and 1% by volume at 4.5$°\mathrm{C}$ [45]. The constraint related to ${\mathrm{H}}_2\mathrm{S}$ content in Cases B and C also influenced ${\mathrm{C}\mathrm{H}}_4$ production, which increased with respect to ${\mathrm{C}\mathrm{H}}_4$ production in Case A. By further adding constraints related to the contents of ${\mathrm{C}\mathrm{H}}_4$, ${\mathrm{H}}_2$, and ${\mathrm{N}\mathrm{H}}_3$ in the produced biogas (Case C), ${\mathrm{C}\mathrm{H}}_4$ and biogas production increased by 12.3% and 1.1%, while the production of ${\mathrm{H}}_2$, ${\mathrm{H}}_2\mathrm{S}$, and ${\mathrm{N}\mathrm{H}}_3$ was reduced by 28.5%, 19.2%, and 80.9% with respect to the initial state.

The variations in biogas compounds for optimized states (Cases A, B, and C) with respect to the initial state are presented in Figure 14. It is evident that in Case C, the content of desired ${\mathrm{C}\mathrm{H}}_4$ increased by 11.1%, while the contents of undesired ${\mathrm{H}}_2$, ${\mathrm{H}}_2\mathrm{S}$, and ${\mathrm{N}\mathrm{H}}_3$ reduced by 30%, 20%, and 81% when compared to the initial state.

The presented results show that temperature, biological additives (bacteria, enzymes), and inorganic additives (which impact the $\mathrm{p}\mathrm{H}$ value), as well as their interaction, play an important role in the improvement of the AD process. Furthermore, the analysis of the obtained results shows that the optimization of the AD process in Case A delivers maximal biogas production, while the optimization in Case C maximizes the quality of the produced biogas.

4. Conclusions

This work investigated the suitability and efficiency of coupling a complex BioModel and an approximation-type gradient-based optimization method into a numerical optimization procedure during the AD process. Numerical experience has shown that the proposed approach provides a reliable and stable numerical procedure capable of optimizing the AD process quite efficiently. Moreover, no problems could be observed by engaging a gradient-based optimizer in combination with design derivatives and obtained by finite differences.

The best quality of biogas produced in reasonably large volumes at an acceptable cost can be reached by simultaneously maximizing biogas production and minimizing the heating energy needed. In this context, one must constrain the $\mathrm{p}\mathrm{H}$ value, the maximal total initial concentrations of bacteria, minimal allowed content of ${\mathrm{C}\mathrm{H}}_4$ and maximal allowed contents of ${\mathrm{H}}_2$, ${\mathrm{H}}_2\mathrm{S}$, and ${\mathrm{N}\mathrm{H}}_3$ in the produced biogas. In such an optimization case, in this work, the content of ${\mathrm{C}\mathrm{H}}_4$ in biogas increased by 11%, while the contents of ${\mathrm{H}}_2$, ${\mathrm{H}}_2\mathrm{S}$, and ${\mathrm{N}\mathrm{H}}_3$ reduced by 30%, 20%, and 81%, respectively; in the optimized system, the bioreactor temperature was $25 °\mathrm{C}$ and the pH value was 7.7.

Since the proposed optimization procedure has shown itself to be numerically stable and efficient, future work should be focused in two directions. Firstly, the AD model should be further improved by considering more precise and various bacterial growth and inhibition models. Secondly, it would be beneficial to widen the AD process’s optimization possibilities by including parameters related to the feeding substrate and feeding strategy.