Full-Scale Digesters: Model Predictive Control with Online Kinetic Parameter Identification Strategy

Citation: Cortés, L.G.; Barbancho, J.; Larios, D.F.; Marin-Batista, J.D.; Mohedano, A.F.; Portilla, C.; de la Rubia, M.A. Full-Scale Digesters: Model Predictive Control with Online Kinetic Parameter Identification Strategy. Energies 2022, 15, 8594. https://doi.org/10.3390/ en15228594 Academic Editor: Giovanni Esposito Received: 21 August 2022 Accepted: 10 November 2022 Published: 16 November 2022 Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations. Copyright: © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/). energies Article Full-Scale Digesters: Model Predictive Control with Online Kinetic Parameter Identification Strategy Luis G. Cortés 1,*, J. Barbancho 1, D. F. Larios 1, J. D. Marin-Batista 2,3, A. F. Mohedano 3, C. Portilla 4 and M. A. de la Rubia 3 1Departamento de Tecnología Electrónica, Escuela Politécnica, Universidad de Sevilla, 41011 Seville, Spain 2Efuels Technologies Ltd., 42-44 Bishopgate, London EC2N 4AH, UK 3Departamento de Ingeniería Química, Campus de Cantoblanco, Universidad Autonoma de Madrid, 28049 Madrid, Spain 4Facultad de Minas, Universidad Nacional de Colombia, Robledo, Medellín 050034, Colombia *Correspondence: luicor[email protected] Abstract: This work presents a nonlinear model predictive control scheme with a novel structure of observers aiming to create a methodology that allows feasible implementations in industrial anaerobic reactors. In this way, a new step-by-step procedure scheme has been proposed and tested by solving two specific drawbacks reported in the literature responsible for the inefficiencies of those systems in real environments. Firstly, the implementation of control structures based on modeling depends on microorganisms’ concentration measurements; the technology that achieves this is not cost-effective nor viable. Secondly, the reaction rates cannot be considered static because, in the extended anaerobic digestion model (EAM2), the large fluctuation of parameters is unavoidable. To face these two drawbacks, the concentration of acidogens and methanogens, and the values of the two reaction rates considered have been estimated by a structure of two observers using data collected by sensors. After 90 days of operation, the error in convergence was lower than 5% for both observers. Four model predictive controller (MPC) configurations are used to test all the previous information trying to maximize the volume of methane and demonstrate a satisfactory operation in a wide range of scenarios. The results demonstrate an increase in efficiency, ranging from 17.4% to 24.4%, using as a reference an open loop configuration. Finally, the operational robustness of the MPC is compared with simulations performed by traditional alternatives used in industry, the proportional-integral- derivative (PID) controllers, where some simple operational scenarios to manage for an MPC are longer sufficient to disrupt a normal operation in a PID controller. For this controller, the simulation shows an error close to the 100% of the reference value. Keywords: anaerobic digestion; asymptotic observer; homogeneous reaction systems; kinetic parameter observer; model predictive control; step-ahead 1. Introduction The main attribute of industrial processes and manufacturing is the large-scale transformation of mass and energy into high-end valuable specific materials and compounds [ 1 ]. In essence, chemical and biological processes are imperative for those purposes because it is primarily at this stage where the elements in a reaction can be manipulated to achieve the desired attributes [ 2 – 4 ]. All the above is possible if a control and monitoring scheme is added to guarantee safe and stable operation [ 5 ]. Specifically, in anaerobic systems, the transformation of substrates into biogas and digestate has attractive advantages in comparison to classical alternatives such as aerobic digestion or composting [ 6 , 7 ]; it returns little sludges, has a positive overall energy balance, and also has enormous potential to reduce challenging and concentrated substrates such as animal wastes [ 8 , 9 ]. However, this is an elaborated biochemical process that includes diverse microorganisms, making the Energies 2022,15, 8594. https://doi.org/10.3390/en15228594 https://www.mdpi.com/journal/energies Energies 2022,15, 8594 2 of 23 process nonlinear [ 10 ]. In this area, it is great of interest to derive models that would be as sensitive as possible to the lack of phenomenological knowledge [ 11 , 12 ]. Thus, an interesting alternative to represent this phenomenon is the extended anaerobic digestion mathematical model 2 (EAM2) [ 13 ], designed for control and monitoring purposes. It circumvents the difficulties in locating the biological lack of knowledge in dedicated terms, namely reaction rates [ 14 ]. Additionally, compared with its predecessor, in the anaerobic digestion mathematical model 2 (AM2), two additionally yield parameters were included to consider the influence of proteins and aminoacids on the dynamic of alkalinity and consider more substrates to be modeled. In order to allocate a strategy that better fits the challenge mentioned above, several techniques have been proposed from different control areas. Some of them are nonlinear proportional-integral-derivative (PID) controllers, linear quadratic regulation (LQR) for nonlinear systems [ 15 ], nonlinear model predictive control (NMPC) [ 16 ], and receding horizon controllers [ 17 , 18 ], to name the most relevant cases reported in the literature. Although these strategies have proven to be effective, the remaining inefficiencies surge from the fact that some of the parameters are not available to be measured online, and the concentrations of the microorganisms cannot be measured periodically because on-site analytical procedures are challenging to perform [ 5 ]. A nonlinear model-based control scheme has been proposed to achieve a satisfactory operation. At the same time, the representation of anaerobic digestion (AD) processes is reduced to a minimal number of assumptions in the model-building exercise, maintaining the performance [ 18 , 19 ]. Among the alternatives registered in the literature, the Model Predictive Control (MPC) [ 17 , 18 ] offers several benefits compared with traditional methods [ 14 , 20 , 21 ]. This algorithm uses a mathematical model to predict and follow the system’s evolution in advance. One of the main advantages is the explicit use of physical and operational constraints in the controller, allowing it to stand around working under feasible operational conditions [22,23]. To achieve this objective, a necessary condition must be satisfied to unlock the door and enable the operation of control schemes based on modeling. All the information needed from the anaerobic reactor must be available periodically, so a group of reliable sensors should be considered to measure, in real-time, the evolution of reactions in the homogeneous reaction system [ 24 ]. Some key variables, the concentrations of microorganisms, as well as the growth rates of reactions (considered variable due to the large fluctuation of parameters in AM2 [ 10 ]), are usually tested using analytical methods through specialized laboratories on-site. The cost and duration of these analyses limit the frequency of the sampling [ 25 ]. Additionally, the available online sensors on the market for biomass or metabolites need to be adequate for regular industrial applications [ 14 ]. The previous assumption sets a significant challenge because collecting this information is critical for the operation of the controller algorithm. In general, the unavailability of specialized instruments incentive the design and construction of observers dedicated to the reaction process. Hence, the design of software sensors, based on the available measurements [ 14 , 26 ], allows to guarantee a streamline of information between the reactor and the control algorithm [ 27 ]. State and parameter observers established over accurate and strict models are the best alternatives to find a way to obtain the lack of the data needed [ 28 – 30 ]. For anaerobic reaction processes, the nonlinear AM2 has been widely proven for monitoring and control purposes [ 10 , 13 ]. This model, constructed under a mass-balance framework, can adapt its dynamics to various AD configurations of reactors and elude the drawback related to the absence of phenomenological knowledge, concentrating this information in non-sensitiveness-dedicated terms named reaction rates [ 10 , 23 ]. The complexity of the biochemical process considers different types of microorganisms, making the process nonlinear and unpredictable, impractical for the traditional online estimation strategies, especially if the operation of the homogeneous system is conditioned by constraints [ 5 ]. To deal with these inconveniences, a sequential structure in cascade has been proposed [ 14 ]. The main idea is to estimate the concentration of acidogens and methanogens in the absence of reaction rates. Then, the next step is to use this information to perform a new observation Energies 2022,15, 8594 3 of 23 that estimates the reaction rates. Acquiring in-depth knowledge of the evolution of reaction rates represents a challenge since it reveals the degree of completion of both reactions during the process [5,31]. This paper is organized as follows. Section 2explains the extended AM2 mathematical model used for control purposes, where additional terms are included due to the influence of proteins and aminoacids on the dynamics of alkalinity. The modification of the AM2 model [ 10 ] unlocks the door to add a broad spectrum of substrates allowed to be used by the control algorithm. Then, Section 3presents the parameter identification procedure based on optimization that uses a novel step-ahead algorithm to considerably improve the efficiency of adjusting the dynamics of the model to the data. The asymptotic state observer that operates without information on the reaction rates is presented to estimate the concentrations of acidogens and methanogens. Next, using the previous information, a second observer estimates the evolution of the two reaction rates considered. In Section 4, the information taken from observers and the reactor’s measurements is used to enable the operation of the MPC. Finally, Section 5shows the improvement achieved by using the MPC controller, measuring the increase in the volume of methane produced. The robustness of the MPC proposed is compared with the PID controller. 2. Optimal Parameter Identification Algorithm The model EAM2, initially proposed by Bernard et al. [ 10 ] as AM2 in 2002 and modified by Cortés et al. [ 13 ] in 2022, was selected to support the new monitoring and control structure presented in this paper. The new model solves two issues documented in the literature as the main drawbacks in implementing automated decision support systems on anaerobic reactors [ 11 ]. First, the necessity to characterize the inlet substrate [ 10 ]. EAM2 considers a wider range of substrates than AM2 because it takes into account the effect of ammonium produced from proteins and aminoacids that affect the dynamics of alkalinity. Second, the concentrations of acidogens and methanogens are estimated using a specific state observer. The above-mentioned solutions allow balancing the absence of reliable sensors or analytical methods to obtain the needed information periodically [10]. 2.1. Proposed Mass-Balance Mathematical Model for Anaerobic Reactors The extended mathematical model used in this paper is based on the AD represented by two main groups of microorganisms. Two consecutive process describe the AD phenomenon on reactors, the acidogenic bacteria (represented by X 1 ) that consume the organic substrate (represented by S 1 ) and produce mainly volatile fatty acids (mainly acid acetic, represented by S 2 ). Thus, the second group of microorganisms, the methanogenic archaea (represented by X2), consumes S2and produces CO2and CH4. Two yield parameters were added to represent the generation of ammonium due to the fermentation of aminoacids or the consumption by the growth of microorganisms. The terms added are K Z,1 and K Z,2 , respectively. The model presented in Cortés et al. [ 13 ] considers the following dynamical equations: dX1 dt =X1(µ1−αD)(1) dX2 dt =X2(µ2−αD)(2) dS1 dt =D(S1in −S1)−k1ψ1S1 KS1+S1(3) dS2 dt =D(S2in −S2) + k2ψ1S1 KS1+S1−k3ψ2S2 KS2+S2(4) Energies 2022,15, 8594 4 of 23 dZ0 dt =D(Zin −Z) + kZ,1ψ1S1 KS1+S1+kZ,2ψ2S2 KS2+S2(5) dC dt =D(Cin −C)−qC+k4µ1X1+k5µ2X2(6) where X 1 is the concentration of acidogens (kg/m 3 ), X 2 is the concentration of methanogens (kg/m 3 ), S 1 is the concentration of organic substrate characterized by chemical oxygen demand (COD) (kg/m 3 ), S 2 is the concentration of volatile fatty acids (VFA) (g acetic acid/L), Zis the total alkalinity (kg/m 3 ), and Cis the concentration of total inorganic carbon (kg/m 3 ). The variables S 1in ,S 2in ,C in , and Zin are, respectively, the influent concentrations of S 1 ,S 2 ,C, and Z. The variable α is related to the reactor design, with values ranging between 0 and 1. The value 0 corresponds to an ideal fixed-bed reactor, and α= 1, to an ideal continuously stirred tank reactor (CSTR). Dis the dilution rate (d −1 ) and is the inverse to the hydraulic retention time, HRT = 1/D, ψ1 is the maximum acidogenic bacteria growing rate (d −1 ), ψ2 is the maximum methanogenic archaea growing rate (d −1 ), K S1 is the half-saturation constant of S 1 (kg/m 3 ), and K S2 is the half-saturation constant of S 2 (g acetic acid/L). The yield parameter k 1 is a constant for substrate degradation, k 2 is a constant for VFA production (mmol/kg), k 3 is a constant for VFA consumption (mmol/kg), k 4 is a constant for CO 2 production (mmol/kg), k5 is a constant for CO 2 production (mmol/kg), and k 6 is a constant for CH 4 production (mmol/kg). The Monod-type Equations (7) and (8) are characterized by the following two reaction rates: µ1=ψ1S1 KS1+S1(7) and, µ2=ψ2S2 KS2+S2(8) BothMonod-typekineticsdescribethegrowth ofacidogenic bacteria ψ1 andmethanogenic archaea ψ2 because, during the fermentation process, the biomass does not register possible VFA accumulation and consequently inhibition. Finally, the methane flow rate produced, qM, is proportional to the reaction rate of methanogenesis, as shown in Equation (9): qM=k6ψ2S2 KS2+S2(9) 2.2. Experimental Data Results The data considered for modeling purposes [ 32 ] were used previously in Cortés et al. [ 13 ] to characterize a CSTR pilot thermophilic digester (150 L) using the AM2 extended mathematical model. The total extension of the experiment was 338 days, where the effect of solid retention time (SRT) and the dynamics of the system in a wide range of operational conditions were studied. However, only a specific range of data was used, where the acidogens and methanogens operated under stable and homogeneous conditions discarding unstable scenarios unfavorable for modeling purposes (see the conditions argued to discard the data). Finally, the data considered represent a total of 207 days. The values of COD at the inlet and VFA and COD at the effluent are shown in Cortés et al. [13]. 3. Optimal Parameter Identification and Online Measurements 3.1. Parameter Identification Using Pattern Search by Step-Ahead Every sampling time, the parameter identification algorithm step-ahead uses the EAM2 to calculate the following value of the system in advance. The procedure established that on the following step, the predicted value is compared with the same variable measured. The difference in magnitude of both represents the cumulative error. The procedure is performed during all the experiment. The following vector Estores all this information. Energies 2022,15, 8594 5 of 23 E= [error1,error2, . . . , errorn](10) nrepresents the total number of days evaluated. Finally, an optimization problem is solved aiming to minimize all the errors calculated along the simulation. min p(k)E subject to: nm(k+1) = f(nm(k),u(k)), (11) nnm(k+1) = f(nnm(k),u(k)), 0⩽p(k)⩽pmax,∀k=1, . . . , tf nm represents all the states values measured S 1 ,S 2 ,C, and Z. nnm represents all the non-measurable states X 1 and X 2 . The parameters considered to be identified are p(k) e{µ1max , µ2max , KS1 , KS2 , k1 , k2 , k3 , k4 , k5 , k6 , KZ1 , KZ2 , Zin} . pmax represents the maximum value affordable by each parameter. This algorithm has already been used in a previous work with the same elements, real data, and phenomenological models [13]. 3.2. Asymptotic Estimator When Reaction Rates Are Unknown The absence of information due to the lack of reliable sensors or analysis results opens up an opportunity to substitute the uncertainty by the design of software sensors. In anaerobic reactors, this technology estimates the concentration of acidogens and methanogens without the information of the reaction rates [ 11 ]. It results in a particular category of observer-denominated asymptotics, because two conditions support its operation; the system is still not exponentially observable, and the reaction kinetics are unknown [ 14 ]. The detailed information of the design is presented in Cortés et al. [13]. In the following, we present the conditions that affect the design of the algorithm: the information of the matrix φ is unknown, the yield coefficients from K are fully known, and the number of q state , the number of measured state variables, is the same or higher than the rank of the matrix K(that is, q state = dim( ξ1)≥ rank(K)). Hence, based on the previous information, the general equation of homogeneous reaction systems described by a general nonlinear state space model is presented: dξ dt =Kφ(ξ,t)−Dξ−Q(ξ) + F(12) where dim( ξ ) = dim(F) = dim(Q) = N, dim( φ ) = M and dim(K) = N × M. Thus, the general nonlinear model Equation (12) can be divided as dξa dt =Kaφ(ξa,ξb)−Dξa−Qa+Fa(13) dξb dt =Kbφ(ξa,ξb)−Dξb−Qb+Fb(14) where the rank of Kis p. The submatrix Ka results from a section of Kwith p × M. The submatrix Kb has the remaining information of K. Finally, the matrices (ξa , ξb) , ( Qa , Qb ), and (F a ,F b) are the corresponding parts of ξ ,Q, and Fcaused by the influence of K a and K b . The previous formulation has the following feature: Zob =A0ξa+ξb(15) There exists a transformation that considers Zob as a linear combination of ξa and ξb as seen in Equation (15). Finally, after processing Equations (13) – (15) , the new state space model is equivalent to: Energies 2022,15, 8594 6 of 23 dξa dt =Kaφ(ξa,ξb)−Dξa−Qa+Fa(16) dZob dt =−DZob +A0(Fa−Qa) + (Fb−Qb)(17) The expression Fa−Qa = 0 means the partition made by Equations (13) and (14) are appropriate since the new dynamics on Zob are independent from K and φ . In Equation (16) , the independence of ˙ ξa from φ is stated. To extend for detailed information, check the complete demonstration in Cortés et al. [13]. Observer Design Using the nonlinear EAM2, the following equations describe the decoupled subsystem conducted by the state variables X1,X2,S1, and S2: ξ=    X1 X2 S1 S2    ,F=    0 0 DS1in DS2in    ,Q=    0 0 0 0    ,K=    1 0 0 1 −k10 k2−k3    ,φ=µ1X1 µ2X2(18) Considering the previous subsystem, the subsequent state equation is structured as follows: d dt    X1 X2 S1 S2    =    1 0 0 1 −k10 k2−k3    φ1 φ2−D    X1 X2 S1 S2    +    0 0 DS1in DS2in    (19) where ξa and ξb group the measurable and non-measurable states, respectively. Thus, using the equations described in Cortés et al. [13], the expression of Zob is as follows: dZob dt =−DZob +A0(Fa −Qa) + (Fb−Qb)(20) where the matrices related to the equation are: ξb=X1 X2,Zob =Zob1 Zob2,A0="1 k10 k2 k1k3 1 k3#,ξa=S1 S2(21) Fa=DS1in DS2in,Fb=0 0(22) with: Zob =Zob1 Zob2=A0ξa+ξb="1 k10 k2 k1k3 1 k3#S1 S2+X1 X2(23) 3.3. Kinetic Parameter Reaction Estimator This section presents an additional tool that complements the lack of measurements in the reactor. The goal is to estimate the kinetic reactions coming from the reaction system. Figure 1shows the structure proposed to obtain all the information needed to feed the mathematical model contained in the control scheme strategy. Energies 2022,15, 8594 7 of 23 Figure 1. The controller and the observer structure proposed for the anaerobic reactor. The discontinuous line represents the information generated by the observers. At the end, the measured variables nm and the estimated ones converge to feed the controller algorithm. In order to design the kinetic estimator algorithm, the following nonlinear equation that represents the system must be considered: ˙ ξ=Kφ−Dφ−Q+F(24) Equation (24) assumes that the coefficients from K are known, while the dilution rate D , the input flows F , and the gas output flows Q are measured in real-time. Additionally, it is assumed that ξ is fully-known because the non-measured dynamics are reconstructed by the previous asymptotic observer. The vector φ is partially known and is divided as follows: φ=Hρ(25) The matrix H contains the information from the known kinetic reactions, and ρ contains the remaining information, the unknown kinetic reactions. Thus, using Equation (25) in Equation (24) results in Equation (26). ˙ ξ=KHρ−Dφ−Q+F(26) The estimation of reaction kinetics, r 1 and r 2 , is equivalent to estimating the entire vector ρ . Hence, following the structure proposed by Bastin et al. [ 14 ], a new dynamic of the system is presented as follows: ˙ ˆ ξ=KH ˆ ρ−Dξ−Q+F−Ω(ξ−ˆ ξ), (27) ˙ ˆ ρ=(KH)TΓ(ξ−ˆ ξ)(28) where ˆ ρ represents the real-time estimation of ρ . When ξ−ˆ ξ∼ = 0 occurs in Equation (28) , it means that ˙ ˆ ρ∼ = 0 results in a perfect tendency process of convergence to desirable reaction kinetics values ˆ ρ∼ =ρ . There is a mandatory condition. Equation (27) that represents the kinetic reaction estimator has an equivalent structure to the Luenberger observer used in homogeneous reaction systems [14]. ˙ ˆ ξ=Aˆ ξ+Bu | {z } 1∗ +Ω(ξ−ˆ ξ)(29) Energies 2022,15, 8594 8 of 23 The underlined term 1 ∗ in Equation (29) ,A ˆ ξ +Bu, is equivalent to the term KH ˆ ρ− D ξ− Q +Fin Equation (27). Thus, the above-mentioned comparison results in Equation (30): ˙ ˆ ξ=KH ˆ ρ−Dξ−Q+F | {z } 1∗ −Ω(ξ−ˆ ξ)(30) In this way, Equation (30) represents the nonlinear system, while Equation (29) is the equivalent linearized system. Observer Design To design the kinetic reaction estimator, only the measured and estimated states on ˙ ξ are needed, while the term ξ−ˆ ξ is used as a reference to monitor and follow the estimator’s performance. ˆ ξ should converge to ξ as soon as possible ( ξ≈ˆ ξ ). The dynamics of ˆ ˙ ξ depends on the evolution of the dynamics of the error ˙ e = (A −Ω L)e. Thus, the references (the known values) that assist in tuning the kinetic reaction estimator are ξ , ˆ ξ , K , and H . For better understanding, the following definitions are proposed: the term e= ξ−ˆ ξ means the observation error, and the term ˜ ρ=ρ−ˆ ρ means the tracking error. Hence, using Equations (31) and (32): ˙ ˆ ξ=KH ˆ ρ−Dξ−Q+F−Ω(ξ−ˆ ξ)(31) ˙ ξ=Kφ−Dξ−Q+F(32) Deriving on both sides of e and ˜ ρ , it results in ˙ e=˙ ξ−˙ ˆ ξ and ˙ ˜ ρ=˙ ρ−˙ ˆ ρ . Taking Equations (31) and (32) and replacing them by ˙ eresults in: ˙ e=Kφ−Dξ−Q+F−KH ˆ ρ+Dξ+Q−F+Ω(ξ−ˆ ξ)(33) ˙ e=Kφ−KH ˆ ρ+Ωe(34) ˙ e=KHρ−KH ˆ ρ+Ωe(35) ˙ e=KH(ρ−ˆ ρ) + Ωe(36) ˙ e=KH ˜ ρ+Ωe(37) In the same way, the dynamic of the tracking error is as follows: ˙ ˜ ρ=−(KH)TΓe+dρ dt (38) Merging Equations (37) and (38) results in the following dynamic system: ˙ e ˙ ˜ ρ="KH ˜ ρ+Ωe −(KH)TΓe+dρ dt #(39) Thus, after organizing the previous equation, it results in the following system: ˙ e ˙ ˜ ρ=ΩKH −HTKTΓ0e ˜ ρ+"0 dρ dt #(40) Replacing the matrices: A=ΩKH −HTKTΓ0(41) V="0 dρ dt #(42) This results in: ˙ e ˙ ˜ ρ=Ae ˜ ρ+V(43) Energies 2022,15, 8594 9 of 23 According to the information described, the estimation procedure is performed as follows: Ka=−k10 k2−k3,Kb=1 0 0 1,Fa=DS1in DS2in,Fb=0 0(44) The next step is to calculate the matrix A0; thus, it means that: A0=−KbK−1 a(45) In the following equations, the system is divided into the dynamic states to be estimated, using the data from the asymptotic observer, Z 1 and Z 2 , and the measured data S 1 and S2. Solving ξbfrom Equation (15), we obtain: ξb=Zob −A0ξa(46) Equation (46) allows to calculate the following matrices: K=Kb Ka,H=1 0 0 1,ξ=  ξb S1 S2 ,F=Fb Fa(47) Finally, all the information needed to supply Equations (27) and (28) are obtained above. 4. Nonlinear Model Predictive Controller (Mpc) As it is shown in Figure 1, the observer structure calculates the non-measurable dynamic states, the concentrations of X 1 and X 2 (with the asymptotic observer), and the kinetic reaction rates r 1 and r 2 (with the reaction rate observer). Once the information from measurements and virtual sensors has been delivered to the controller, the algorithm is ready to start calculating the control actions after receiving the rest of the information: the operational and physical instructions, as well as the control objectives. 4.1. Controller Design The subsequent conditions emerge as the required rules to take into consideration to design the controller. • The reactor has to be balanced to react against perturbations, always working within the physical and operational boundaries. • Methane production needs to be maximized all the time. • The environmental regulations and the capacity of the reactor to reduce the concentration of substrates at the inlet, S 1 and S 2 , conditioned the programming. Therefore, the rule S 1 (t) + S 2 (t) ≤ Kt d , where Kt d denotes the maximum effluent concentration of the both substrates considered, has to be followed. • The reactor has to be protected against failures due to unexpected high variations on VFA, and in consequence, inhibitions on the metabolism of microorganisms. Thus, an alternative is to use alkalinity as a base to neutralize the level of acids. S2(t) Z(t)=λ(48) Equation (48) tries to preserve a relation between VFA and Zduring the operation of the reactor. According to the literature, λ has to be within 0.1 and 0.3, to maintain the microorganisms in a zone of comfort in regard to the level of acids [11]. 4.2. Structure of the Controller Finally, the structure of the controller is shown in Figure 2. The dilution rate D controls the amount of substrate at the inlet. It results in u C , the variables used by the MPC to manipulate the dynamics of the reactor. u i represents the non-controller variables considered as inlet: S 1in ,S 2in ,C in , and Z in .u pert holds the information about the inputs considered as disturbances. y contains the information of the variables available to be measured. n m are the state dynamics that can be measured: S 1 ,S 2 ,Z, and C. The value of pH Energies 2022,15, 8594 16 of 23 5.4. Mpc Controller without Restrictions Figure 7shows the results of maximizing the methane produced by the reactor using an MPC without the operational and physical restrictions shown in Equation (49). (a) (b) (c) (d) (e) (f) Figure 7. MPC controller without restrictions, maximizing the volume of CH 4 . Time course of ( a ) volume of CH 4 ; ( b ) dilution rate D; ( c ) value of pH; ( d ) dynamic of VFA; ( e ) values of λ ; ( f ) dynamics of S2and Z. Figure 7a shows a substantial increase in the volume of methane produced (see discontinuous red line), where the black line acts as the reference (volume of methane produced due to profile inputs used in parametric optimization), to measure the improvements on efficiency. The correspondent control actions calculated are shown in Figure 7b. The absence of restrictions in the algorithm allows the controller to calculate control actions, neglecting strong changes in operation that could affect the stability of the reactor. As can be seen, maximizing the objective turns the values of D extreme; showing oscillations from one Energies 2022,15, 8594 17 of 23 extreme to the other. Several times the value of D decreases to zero and suddenly rise to the maximum. In the same way, pH had severe changes even beyond the operational limits (see Figure 7c). The relationship between pH and VFA is clear in Figure 7c,d showing a correlation along the simulation. This algorithm exposes high external risks in operation because there exists a lack of information over the system’s operation. Figure 7e shows the corresponding values of λ along the simulation; it appears that the values go beyond the limits repeatedly. As it seems in Figure 7f, the value of alkalinity Z exceeds the value of S2 where at this point, the value of λ goes beyond 1, meaning the reaction system runs inadequately. Given the previous results in Figure 7, the maximization of methane without restrictions leads the system to operate beyond non-feasible limits. Therefore, it is necessary to include extra information about the process to avoid inhibitors along the simulation. Figure 8shows the results obtained due to the operation of the MPC without the restriction of buffering capacity ( λ ) and without the multistart function. Figure 8a shows the volume of methane produced by both the reference (black line), and the MPC (red discontinuous line). A significant increase in the volume of methane produced compared to the reference is observed. Figure 8b shows the control actions calculated by the MPC. As it is shown, the oscillation is severe, however, the results are reliable. Figure 8c,d shows the results of the pH and S2 . As expected, the oscillation of values is intense around the reference (black line) but not attractive for control purposes. Finally, Figure 8e shows the evolution of the non-restricted parameter λ . As expected, due to the absence of restrictions on programming, the maximum value of λhas been exceeded once. (a) (b) (c) (d) Figure 8. Cont. Energies 2022,15, 8594 18 of 23 (e) (f) Figure 8. MPC controller without restrictions. ( a ) Volume of CH 4 ; ( b ) dilution rate D; ( c ) value of pH; (d) dynamic of VFA; (e) values of λ; (f) dynamics of S2and Z. 5.5. Model Predictive Controller with Restrictions The information about operational and physical restrictions is added to the MPC aiming to maintain the microorganisms under feasible operational conditions. As shown in Figure 9a, the discontinuous red line shows the volume of methane produced. The efficiency is maintained despite the restrictions; therefore, this option better fits the system under stable conditions. The correspondent control actions calculated by the optimizer (see Figure 9b) shows a slight attenuation compared with previous simulations. The results of pH in Figure 9c show that the range of values becomes closer to the reference than the previous results. Compared with the previous results, the acidification level is much lower (see Figure 9d). In Figure 9e, the discontinuous blue line represents the limit stated for a regular operation. As it is shown, the discontinuous red line never surpasses the limits beyond λ = 0.8; in some cases, it gets closer. However, it is because the reactor is producing a high amount of methane closer to the operational limits. Figure 9f shows the evolution of the balance between VFA and Z. It is clear that the buffer capacity is always under control. Although the algorithm is restricted, the amount of methane CH 4 produced is higher than the reference. Finally, Figure 10 shows similar results to those shown in Figure 8but programmed without the multistart function. Despite the decrease in efficiency, the results show an adequate operation. (a) (b) Figure 9. Cont. Energies 2022,15, 8594 19 of 23 (c) (d) (e) (f) Figure 9. MPC controller with restrictions, maximizing the volume of CH 4 . ( a ) Volume of CH 4 ; (b) dilution rate D; (c) value of pH; (d) dynamic of VFA; (e) values of λ; (f) dynamics of S2and Z. (a) (b) (c) (d) Figure 10. Cont. Energies 2022,15, 8594 20 of 23 (e) (f) Figure 10. MPC controller with restrictions. ( a ) CH 4 evolution; ( b ) D values calculated by optimizer; (c) pH values; (d) VFA values; (e)λover the process; (f) evolution of S2and Z. Table 1shows a condensed vision of the previous results. Two same alternatives were simulated without the multistart function. As expected, there is a reduction in efficiency because the algorithm chose a random initial point to start the process. The black line is used as a reference to compare easily the results from the methane produced by the MPC [ 13 ]. In the case of the MPC controller that operates without operational restrictions, the improvement achieved was 17.4% and 24.4%, respectively, for the algorithm with no multistart and the multistart function. The reactor that operates with operational restrictions shows an improvement of 18.8% and 20.8% for both MPC schemes, without and with the multistart function, respectively. Although the performance is reduced significantly because of the insertion of restrictions, the continuous operation of the reactor is not broken and can run permanently. Table 1. Improvement in the methane produced by each MPC control scheme. Schemes Without Restrictions With Restrictions No Multistart Multistart No Multistart Multistart Increase on efficiency 17.4 % 24.4% 18.8% 20.9% 6. Conclusions A structure of two observers demonstrates to be effectively solving two specific drawbacks responsible for the inefficiencies implementing control schemes based on modeling. The Asymptotic observer that operates in the absence of the reaction rates demonstrates an error of convergence less than 5% estimating the concentrations of acidogens and methanogens. Then, the previous information enables an upstream operation feeding the kinetic observer that estimates the reaction rates with an error of convergence lower than 5%. In the following, the rest of the parameters, considered as static due to their correlation with the specific biochemical characteristics of the organic matter, were calculated using algorithms based on optimization. All the above-mentioned aspects represent the key that makes it possible to operate MPCs in the industry. Four MPC configurations were used to test the observer and identification structure proposed. The objective was to maximize the volume of methane. While the controller without the operational restrictions showed a considerable increase in efficiency, after operational limits were added to the algorithm, the efficiency of the controller slightly decreased. However, its reduction on efficiency is affordable because the algorithm helped the system to stay in operation continuously. Without operational restrictions, the efficiency reached 24.4%, while the restricted algorithm showed an efficiency of 20.9%. The exhaustiveness in the methodology proposed is worth using because the added information guarantees robustness in operation and an increase in efficiency compared with traditional control alternatives. Energies 2022,15, 8594 21 of 23 Author Contributions: Conceptualization, L.G.C., J.B., D.F.L., J.D.M.-B., C.P. and M.A.d.l.R.; Data curation, L.G.C., C.P. and M.A.d.l.R.; Formal analysis, D.F.L., J.D.M.-B., C.P. and M.A.d.l.R.; Investigation, L.G.C., J.B., D.F.L., J.D.M.-B., A.F.M., C.P. and M.A.d.l.R.; Methodology, L.G.C., J.B. and C.P.; Project administration, J.B. and C.P.; Resources, J.D.M.-B. and A.F.M.; Software, L.G.C.; Supervision, D.F.L., J.D.M.-B. and M.A.d.l.R.; Validation, D.F.L. and M.A.d.l.R.; Visualization, L.G.C.; Writing— original draft, L.G.C.; Writing—review & editing, J.B., D.F.L. and C.P. All authors have read and agreed to the published version of the manuscript. Funding: The authors wish to express their gratitude to Fundación FIDETIA (G91045419), Universidad de Sevilla, Spain, for funding this work. The author Luis G. Cortés express special acknowledgements to Fundación Centro de Estudios Interdisciplinarios Básicos y Aplicados (CEIBA), Colombia, for funding. Laboratorio de Investigaciones en Catálisis y nuevos materiales LICATUC (Universidad de Cartagena). Data Availability Statement: The initial experimental results developed in this study are in the next link: https://doi.org/10.1016/j.procbio.2005.03.073 (accessed on 20 August 2022). Conflicts of Interest: The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Nomenclature C,Cin Total inorganic carbon concentration (kg/m3) DDilution rate (d−1) k1Yield substrate degradation k2Yield for VFA production k3Yield for VFA consumption k4Yield for CO2production k5Yield for CO2production k6Yield for CH4production (m3) KS1Half-saturation constant (kg/m3) KS2Half-saturation constant (kg/m3) KZ,1 Yield for ammonium production (kg/m3) KZ,2 Yield for ammonium production (kg/m3) qCCarbon dioxide flow rate (m3/d) qMMethane flow rate (m3/d) r1, r2Reaction rates (d−1) S1,S1in Organic substrate concentration (kg acid acetic/m3) S2,S2in Volatile fatty acids concentration (kg/m3) X1Concentration of acidogenic bacteria (kg/m3) X2Concentration of methanogenic bacteria (kg/m3) Z,Zin Total alkalinity (kg/m3) αFraction of bacteria in the liquid phase µ1Specific growth rate of acidogenic bacteria (d−1) µ2Specific growth rate of methanogenic bacteria (d−1) ψ1Maximum acidogenic bacteria growth rate (d−1) ψ2Maximum methanogenic bacteria growth rate (d−1) nmVector of measured state variables nnm Vector of non-measurable state variables p,pmax Vector of parameters u,umin,umax Control action (d−1) ξVector of state variables KMatrix with the kinetics of the biochemical φMicrobiological reactions involved on system QGaseous rate of mass outflow from reactor FMass feed rate due to external substrates ρUnknown functions of all states HKnown functions of all the states ΓGain matrix of the updating law ΩSquare gain convergence matrix Energies 2022,15, 8594 22 of 23 λBuffer capacity NpPrediction horizon (d) NuControl horizon (d) References 1. Singh, R.; Sahay, A.; Karry, K.M.; Muzzio, F.; Ierapetritou, M.; Ramachandran, R. Implementation of an advanced hybrid MPC–PID control system using PAT tools into a direct compaction continuous pharmaceutical tablet manufacturing pilot plant. Int. J. Pharm. 2014,473, 38–54. [CrossRef] [PubMed] 2. le Roux, J.; Olivier, L.; Naidoo, M.; Padhi, R.; Craig, I. Throughput and product quality control for a grinding mill circuit using non-linear MPC. J. Process Control 2016,42, 35–50. [CrossRef] 3. Corbett, B.; Macdonald, B.; Mhaskar, P. Model Predictive Quality Control of Polymethyl Methacrylate. IEEE Trans. Control Syst. Technol. 2015,23, 687–692. [CrossRef] 4. Lora Grando, R.; de Souza Antune, A.M.; da Fonseca, F.V.; Sánchez, A.; Barrena, R.; Font, X. Technology overview of biogas production in anaerobic digestion plants: A European evaluation of research and development. Renew. Sustain. Energy Rev. 2017 , 80, 44–53. [CrossRef] 5. Marquez-Ruiz, A.; Mendez-Blanco, C.; Porru, M.; Özkan, L. State and Parameter Estimation Based On Extent Transformations. Comput. Aided Chem. Eng. 2018,44, 583–588. 6. Gaida, D.; Wolf, C.; Bongards, M. Feed control of anaerobic digestion processes for renewable energy production: A review. Renew. Sustain. Energy Rev. 2017,68, 869–875. [CrossRef] 7. Association, E.B. EBA Statistical Report; European Biogas Association: Brussels, Belgium, 2017. 8. Méndez-Acosta, H.; Palacios-Ruiz, B.; Alcaraz-González, V.; González-Álvarez, V.; García-Sandoval, J. A robust control scheme to improve the stability of anaerobic digestion processes. J. Process. Control 2010,20, 375–383. [CrossRef] 9. Anukam, A.; Mohammadi, A.; Naqvi, M.; Granström, K. A Review of the Chemistry of Anaerobic Digestion: Methods of Accelerating and Optimizing Process Efficiency. Processes 2019,7, 504 . [CrossRef] 10. Bernard, O.; Hadj-Sadok, Z.; Dochain, D.; Genovesi, A.; Steyer, J. Dynamicals Model Development and Parameter Identification for an Anaerobic Wastewater Treatment Process. Biotechnol. Bioeng. 2001,75, 424–438. [CrossRef] 11. Kil, H.; Li, D.; Xi, Y.; Li, J. Model Predictive Control with On-line Model Identification for Anaerobic Digestion Processes. Biochem. Eng. J. 2017,128, 63–75. [CrossRef] 12. López Buriticá, K.; Trujillo, S.; Acosta-Medina, C.; Granada D.H. Dynamical Analysis of a Continuous Stirred-Tank Reactor with the Formation of Biofilms for Wastewater Treatment. Math. Probl. Eng. 2015,2015, 512404. [CrossRef] 13. Cortés, L.G.; Barbancho, J.; Larios, D.F.; Jose Marín-Batista.; Mohedano, A.F.; Portilla, C.; de la Rubia, M.A. Full-Scale Digesters: An Online Model Parameter Identification Strategy. Energies 2022,15, 7685. [CrossRef] 14. Bastin, G.; Dochain, D. On-line Estimation and Adaptive Control of Bioreactors. Process Meas. Control. 1990 ,10, 707–723. [CrossRef] 15. Das, S.; Pan, I.; Halder, K.; Das, S.; Gupta, A. LQR based improved discrete PID controller design via optimum selection of weighting matrices using fractional order integral performance index. Appl. Math. Model. 2013,37, 4253–4268. [CrossRef] 16. Mauky, E.; Weinrich, S.; Jacobi, H.; Naegele, H.J.; Liebetrau, J.; Nelles, M. Model Predictive Control for Demand-Driven Biogas Production in Full Scale. Chem. Eng. Technol. 2016,39, 652–664. [CrossRef] 17. González, A.; Adam, E.; Marchetti, J. Conditions for offset elimination in state space receding horizon controllers: A tutorial analysis. Chem. Eng. Process. Process Intensif. 2008,47, 2184–2194. [CrossRef] 18. Hanema, J.; Lazar, M.; Tóth, R. Tube-based LPV Constant Output Reference Tracking MPC with Error Bound. IFAC-Pap. 2017 , 50, 8612–8617. 19. Wang, Y.; Witarsa, F. Application of Contois, Tessier, and first-order kinetics for modeling and simulation of a composting decomposition process. Bioresour. Technol. 2016,220, 384–393. [CrossRef] 20. Andrews, J.F. A Mathematical Model for the Continuous Culture of Microorganisms Utilizing Inhibitory Substrates. Biotechnol. Bioeng. 1968,10, 707–723. [CrossRef] 21. Rossi, E.; Pecorini, I.; Ferrara, G.; Iannelli, R. Dry Anaerobic Digestion of the Organic Fraction of Municipal Solid Waste: Biogas Production Optimization by Reducing Ammonia Inhibition. Energies 2022,15, 5515. [CrossRef] 22. Haldane, J. Enzynmes; MIT Press: Cambridge, MA, USA, 1965; p. 184. 23. Marquez-Ruiz, A.; Mendez-Blanco, C.; Ozcan, L. Constrained Control and Estimation of Homogeneous Reaction Systems Using Extent-Based Linear Parameter Varying Models. Ind. Eng. Chem. Res. 2020,59, 2242–2251. [CrossRef] 24. García-Diéguez, C.; Bernard, O.; Roca, E. Reducing the Anaerobic Digestion Model No. 1 for its application to an industrial wastewater treatment plant treating winery effluent wastewater. Bioresour. Technol. 2013,132, 244–253. [CrossRef] [PubMed] 25. Chaib Draa, K.; Zemouche, A.; Alma, M.; Voos, H.; Darouach, M. A discrete-time nonlinear state observer for the anaerobic digestion process. Int. J. Robust Nonlinear Control 2019,29, 1279–1301. [CrossRef] 26. Song, Y.J.; Oh, K.S.; Lee, B.; Pak, D.W.; Cha, J.H.; Park, J.G. Characteristics of Biogas Production from Organic Wastes Mixed at Optimal Ratios in an Anaerobic Co-Digestion Reactor. Energies 2021,14, 6812. [CrossRef] 27. Batstone, D.; Keller, J.; Newell, B.; Newland, M. Model Development and Full Scale Validation for Anaerobic Treatment of Protein and Tat Based Wastewater. Water Sci. Technol. 1997,36, 423–431. [CrossRef] Energies 2022,15, 8594 23 of 23 28. Isaza-Hurtado, J.; Botero-Castro, H.; Alvarez, H. Robust Estimation for LPV Systems in the Presence of Non-uniform Measurements. Automatica 2020,115, 108901. [CrossRef] 29. Meegoda, J.N.; Li, B.; Patel, K.; Wang, L.B. A Review of the Processes, Parameters, and Optimization of Anaerobic Digestion. Int. J. Environ. Res. Public Health 2018,15, 2224. [CrossRef] 30. Bora, B.J.; Dai Tran, T.; Prasad Shadangi, K.; Sharma, P.; Said, Z.; Kalita, P.; Buradi, A.; Nhanh Nguyen, V.; Niyas, H.; Tuan Pham, M.; et al. Improving combustion and emission characteristics of a biogas/biodiesel-powered dual-fuel diesel engine through trade-off analysis of operation parameters using response surface methodology. Sustain. Energy Technol. Assess. 2022 ,53, 102455. [CrossRef] 31. Condrachi, L.; Vilanova, R.; Meneses, M.; Barbu, M. Anaerobic Digestion Process Control Using a Data-Driven Internal Model Control Method. Energies 2021,14, 6746. [CrossRef] 32. de la Rubia, M.; Perez, M.; Romero, L.; Sales, D. Effect of Solids Retention Time (SRT) on Pilot Scale Anaerobic Thermophilic Sludge Digestion. Process Biochem. 2006,41, 79–86. [CrossRef] 33. Valencia, F.; López, J.D.; Núñez, A.; Portilla, C.; Cortes, L.G.; Espinosa, J.; Schutter, B.D. Congestion Management in Motorways and Urban Networks through a Bargaining-Game-Based Coordination Mechanism; Springer: Berlin/Heidelberg, Germany, 2015; pp. 1–40. 34. Batstone, D.; Keller, J.; Angelidaki, I.; Kalyuzhnyi, S.; Pavlostathis, S.; Rozzi, A.; Sanders, W.; Siegrist, H.; Vavilin, V. Anaerobic Digestion Model No 1 (ADM1). Water Sci. Technol. 2002,45, 65–73. [CrossRef] [PubMed]