- Email: [email protected]
Anaerobic Digestion of Microalgae: Identification for Optimization and Control
Proceedings of the 18th World Congress The International Federation of Automatic Control Milano (Italy) August 28 - September 2, 2011
Anaerobic Digestion of Microalgae: Identification for Optimization and Control ⋆ Elliot Cameron ∗ Francis Mairet ∗∗ Olivier Bernard ∗∗ Monique Ras ∗∗∗ Laurent Lardon ∗∗∗ Jean-Philippe Steyer ∗∗∗ Benoˆıt Chachuat ∗,∗∗∗∗ ∗
Department of Chemical Engineering, McMaster University, Canada (Email: [email protected]) ∗∗ COMORE, INRIA Sophia-Antipolis, France (Email:{francis.mairet,olivier.bernard}@inria.fr) ∗∗∗ LBE, INRA, Narbonne, France (Email:{rasm,lardonl,jean-philippe.steyer}@supagro.inra.fr) ∗∗∗∗ Department of Chemical Engineering, Imperial College London, UK (Email: [email protected]) Abstract: Coupling an anaerobic digester to a microalgal culture is currently considered one of the most promising avenues towards the production of renewable bioenergy, either in the form of biodiesel or biogas. Accurate mathematical models are crucial tools to assess the potential of such coupled biotechnological processes and help optimize their design, operation and control. This paper focuses on the compartment of anaerobic digestion of microalgae. Using experimental data for the anaerobic digestion of Chlorella vulgaris, a grey-box model is developed that allows good prediction capabilities and retains low complexity. The proposed methodology proceeds in two steps, namely a structural and a parametric identification steps. The fitted model is then used to conduct preliminary optimization for the production of biogas from Chlorella vulgaris. The results provide some insight into the potential for bioenergy production from the digestion of microalgae and, more generally, the coupled process. Keywords: Bioprocess, anaerobic digestion, microalgae, model identification, optimization 1. INTRODUCTION Microalgae are currently considered one of the most promising feedstocks for biofuels (Chisti, 2008; Wijffels and Barbosa, 2010). But on the path to making large-scale microalgae culture sustainable, one needs to consider the management of large quantities of residual biomass along with the supply of large amounts of fertilizers. Not only does anaerobic digestion appear to be in an ideal position for addressing those challenges, but it also presents very favorable economic and energetic performance (Sialve et al., 2009; Mussgnug et al., 2011). The anaerobic digestion ecosystem is notoriously complex, involving hundreds of bacterial species. While many models of the anaerobic digestion processes have been developed since the 1970s, including the well-accepted ADM1 model (Batstone et al., 2002), there has been little direct study on the modeling of the anaerobic digestion of algal biomass. In this context, special attention must be paid to the low biodegradability of common microalgae and the large nitrogen content of algal biomass which, when converted into ammonia, may inhibit bacterial activity (Koster and Lettinga, 1984). The focus in this paper is on modeling of the anaerobic digestion of microalgae. Our objective is to develop and identify grey-box models that allow good prediction capabilities, while keeping their complexity as low as possible for use in optimization and control. ⋆ Elliot Cameron is grateful for graduate scholarship from NSERC and OGS. Support from the McMaster Advanced Control Consortium (MACC) is also gratefully acknowledged.
978-3-902661-93-7/11/$20.00 © 2011 IFAC
Modeling biotechnological processes is a tricky task, for which systematic methods and tools are still lacking. In practice, greybox models are often considered, which consist of two parts (Bernard and Queinnec, 2008): a first-principles part derived from mass-balance considerations; and a phenomenological part describing the biological reactions. The model development and identification methodology presented in this paper considers such a model structure and proceeds in two main steps. First, the minimum number of reactions needed to match the variability in a given experimental data set is determined. This allows selection of a consistent reaction scheme, which in turn imposes the key state variables. Next, the mass-balance equations and phenomenological laws are formulated, and the unknown stoichiometric/kinetic parameter values are estimated in order for the model predictions to match the experimental data. In particular, a systematic parametric identification procedure is proposed, whereby minimal parameter combinations are determined on account of both the fit quality and the estimated parameter confidence. This methodology is applied to a set of data obtained over 140 days on an experimental digester fed with Chlorella vulgaris (Ras et al., 2011), and the reaction scheme is based on the recent work by Mairet et al. (2011). An application of the calibrated model is for optimizing the production of biogas from Chlorella vulgaris. One such preliminary optimization study is conducted in this paper, thereby gaining some insight into the potential for bioenergy production from the digestion of microalgae and, by extension, from integrated biotechnological processes that couple a microalgal culture with an anaerobic digester.
5052
10.3182/20110828-6-IT-1002.03196
18th IFAC World Congress (IFAC'11) Milano (Italy) August 28 - September 2, 2011
where ξ, ξ in ∈ IRnξ are the state variables in the reactor liquid phase and in the feed stream; D, the dilution rate; K ∈ IRnξ ×nr , the stoichiometric matrix; r ∈ IRnr , the reaction rate vector; and q ∈ IRnξ , the gaseous rate vector.
−1
Dilution rate (d )
0.4 0.3 0.2 0.1
−1
Substrate additions (gCOD.L )
0
3.1 Determining Model Complexity 0
20
40
60
80
100
120
140
160
The first step consists in identifying the reaction network required to explain the variability of a set of experimental data. Determining the minimum number of reactions required to explain the set of experimental data can be performed according to the principal component analysis (PCA)-based methodology outlined in (Bernard and Bastin, 2005; Bernard et al., 2006).
6
4
2
0
0
20
40
60
80 time (d)
100
120
140
160
Fig. 1. Feed conditions of the anaerobic digester The remainder of the paper is organized as follows. The experimental setup and data are briefly described in §2. Then, the structural and parametric identification steps are detailed in §3 and §4, respectively, and the results obtained for the experimental data are presented. Finally, the fitted model is used to conduct preliminary optimization for the production of biogas from microalgae in §5. 2. EXPERIMENTAL METHODOLOGY Experimental Setup An experimental study of the anaerobic digestion of Chlorella vulgaris was performed over 150 days. This study involved the operation of a 1 L continuous-stirredtank anaerobic digester, with 0.1 L headspace, maintained at 35 ◦ C, and with no pH control. Feed Stream Characterization A detailed characterization of the feed stream can be found in (Ras et al., 2011; Mairet et al., 2011). Feed was introduced as slugs, and equal volumes of reactor medium were removed to maintain a constant liquid volume. A daily average of the dilution rate along with the substrate additions are shown in Fig. 1. The average concentrations of other relevant inlet components are provided in Tab. 1. Table 1. Feed Stream Composition Volatile fatty acids Inorganic nitrogen Inorganic Carbon Inert Cations - Anions pH
0.0 gCOD/L 0.011 mol/L 0.017 mol/L 0.017 mol/L 9.4
Measurement Techniques Readings of the following quantities were taken on an approximately daily basis: total COD (by colorimetric method); ion concentrations (by ion chromatography); biogas volume (by water displacement); biogas composition (by gas chromatography); and pH (by colorimetric method). Random samples were also selected on a less frequent basis for volatile fatty acid (VFA) determination. See (Ras et al., 2011) for more details. 3. STRUCTURAL MODEL IDENTIFICATION A general mass-balance model of the following form is considered for the system: ˙ ξ(t) = Kr (·) + D(t) ξ in (t) − ξ(t) − q (·) , (1)
The number of state variables that can be directly measured is quite limited in most biotechnological processes, and the available measurements often correspond to combinations of the state variables, such as total COD. Sometimes, measurements can also be related to the reaction rates, for example gaseous outflows of low soluble species such as methane. The analysis throughout this subsection is thus conducted for an output vector of the form: ′ C′ ξ y , (2) = y= C′′ r (·) y′′ where ny = n′y + n′′y is the total number of measurements, such ′ ′′ that ny ≥ nr . The matrices C′ ∈ Rny ×nξ and C′′ ∈ Rny ×nr represent the combinations for the states and reaction rates. For any time instants ti < tj , let u(·, ·), v(·, ·) and w(·, ·) be defined as: Z tj ′ in ′ ′ ′ D(t)[y (t) − y (t)]dt y (t ) − y (t ) − j i ti u(ti , tj ) = Z tj ′′ y (t)dt t Z tj Z tj i q(·)dt. r(·)dt, and w(ti , tj ) = v(ti , tj ) = ti
ti
Suppose that N records of u(·, ·), v(·, ·) and w(·, ·) with N > nr are available, and define the matrices U = (u(t0 , t1 ), . . . , u(tN −1 , tN )) and V, W alike. Then, from (1) and (2), the following linear relationship holds: ′ ′ CK C U= V+ W. (3) 0 C′′ | {z } =: Γ
Under the additional assumptions that the matrix Γ is full rank—which requires that the measurements y be nonredundant—and that the measured liquid species y′ are not involved in liquid-gas transfer— which enforces W = 0—the number of reactions required to describe the data can be directly be assessed from the PCA of U. Specifically, each principal axis is representative of a given reaction, and the corresponding principal component represents the relative proportion of the overall variability in the data set that can be accounted for by that reaction. Dividing by the largest eigenvalue makes it easy to determine cumulative variability, and thus the minimum number of reactions required to account for a variability in the data greater than a given threshold; e.g., 95% of the variability. In making this PCA analysis, care should always be taken to normalize and center the data in each row of U as:
5053
18th IFAC World Congress (IFAC'11) Milano (Italy) August 28 - September 2, 2011
µ1 (S1 ) = µ ¯1
100
S1 , S1 + KS1
cumulated variability [%]
90
µ3 (S3 , NH3 ) = µ ¯3
80 70
50 40 30 20 10
0
1
2
3
4
number of reactions
Fig. 2. PCA results e i = U√i − µ(Ui ) , U N σ(Ui )
i = 1, . . . , ny .
Application to Anaerobic Digestion of Microalgae The matrix U is calculated using the available liquid-phase measurements in total COD, inorganic nitrogen and VFA concentrations, along with the methane flow rate measurements that are representative of the methanogenic bacteria activity; moreover, the time instants t1 , t2 , . . . , tN +1 are directly taken as the measurement times, and cubic spline interpolants are used to compute the integral terms. Note that none of the measured liquid-phase species are involved in liquid-gas transfer, and therefore the PCA results can be interpreted in terms of biological reactions. PCA is applied to the centered and normalized data, and the cumulative variability for the principal components is presented in Fig. 2. These results indicate that a minimum of 2 or 3 reactions are sufficient for explaining, respectively, 95% and 99% of the variability in the data. For the purposes of this work a three reaction scheme is considered subsequently. 3.2 Three-Reaction Model Having settled on the number of reactions required to describe the available data, the next step is to determine what three reactions adequately perform this function. It should be noted that much individual expertise associated with grey box modeling may come into play at this stage, and this choice can be somewhat trial-and-error. In this work, the three-reaction model recently developed by Mairet et al. (2011) is considered. This model involves two hydrolysis+acetogenesis steps in parallel, whereby sugars+lipids (S1 ) and protein (S2 ) are converted into VFA (S3 ), followed by a methanogenic step. It is assumed that each reaction is performed by specific bacterial populations X1 , X2 and X3 . • Hydrolysis+acetogenesis of sugars+lipids: µ1 (·) X1
α1 S1 + α2 NH+ −→ X1 + α3 S3 + α4 CO2 4 • Hydrolysis+acetogenesis of proteins: µ2 (·) X2
α5 S2 −→ X2 + α6 S3 + • Methanogenesis: µ3 (·) X3
α7 NH+ 4
+ α8 CO2
α9 S3 + α10 NH+ −→ X3 + α11 CH+ 4 4 + α12 CO2
The biological kinetics are taken as:
S3 S3 + KS3 +
S32 KI3
(4)
Observe that knowledge of the pH is required to compute the dissolved NH3 concentration in µ3 . The following chargebalance equation assumes that all acid/base pairs are in equilibrium with h := 10−pH , and that the pH range of operation is lower than 8: h KC Z +h+ N− C h + KN KC + h KVFA− K H2 O − S3 − = 0, (5) γVFA (KVFA− + h) h with C and N denoting total inorganic carbon and nitrogen, respectively; KC , KN , KVFA and KH2 O , the dissociation con+ − stants for HCO− 3 /CO2 , NH3 /NH4 , VFA /HVFA and water, respectively; Z, the difference between inert cation and anion charges; and γVFA = 64 gCOD /mol, by assimilating VFA to pure acetate.
60
0
S2 , S2 + KS2 KINH3 . KINH3 + NH3
µ2 (S2 ) = µ ¯2
Under perfect mixing conditions, the mass-balance equations for the state variables S1 , X1 , S2 , X2 , S3 , X3 , N , C, and Z, along with the concentration Si of inert COD, read: S˙ 1 = D(β1 S in − S1 ) − α1 µ1 X1 (6) X˙ 1 = (µ1 − D)X1 S˙ 2 = D(β2 S in − S2 ) − α5 µ2 X2 X˙ 2 = (µ2 − D)X2
S˙ 3 = − D S3 + α3 µ1 X1 + α6 µ2 X2 − α9 µ3 X3 X˙ 3 = (µ3 − D)X3 N˙ = D(N in − N ) − α2 µ1 X1 + α7 µ2 X2 − α10 µ3 X3
C˙ = D(C in − C) + α4 µ1 X1 + α8 µ2 X2 + α12 µ3 X3 − ρCO2 Z˙ = D(Z in − Z) S˙ i = D(βi S in − Si ).
This formulation neglects the bacterial decay rates. It also assumes negligible NH3 and VFA gaseous losses. Expressions of the gas outflow rate and composition are given in (Mairet et al., 2011). 4. PARAMETRIC MODEL IDENTIFICATION
Once a candidate model structure has been identified, adequate values must be given to the parameters in the corresponding equation system. While accurate values for some of these parameters can be found in the literature, such as acid/base dissociation constants, other parameters may be more uncertain and need to be estimated from experimental data; this typically includes a large number of stoichiometric and kinetic parameters as well as initial conditions for the state variables. 4.1 Systematic Parameter Identification Procedure The heart of parametric model identification is to determine estimates for the uncertain parameters, which are denoted by p subsequently. Such estimates can be obtained by minimizing a specified objective function, typically a weighted sum of the square errors between measured outputs ym and corresponding model predictions y at a number of observation times tj , j = 0, . . . , N :
5054
18th IFAC World Congress (IFAC'11) Milano (Italy) August 28 - September 2, 2011
min J := p
ny N X X i=1 j=0
2 ωi,j yi (tj ; p) − yjm (tk ) ,
(7)
where ωi,j denotes the weighting coefficient for the jth measurement of the ith output; for example, these weights can be chosen as the inverse of the standard deviation of the corresponding output measurement. Moreover, the model outputs y correspond to some state variables ξ or, more generally, functions of the state variables (see §3.1). Mathematical models for biotechnological processes often comprise many uncertain parameters and relatively few measured outputs, which makes them notoriously hard to calibrate due to structural/practical identifiability issues (Dochain and Vanrolleghem, 2001). For this reason, attempting to fit all the parameters simultaneously almost always results in nearly zero confidence in the estimated parameters. Instead, two typically conflicting objectives must be accommodated in calibrating biotechnological process models, namely the fit quality and the estimated parameter confidence. Obtaining a close fit between the model predictions and the experimental data is always favored by increasing the number of degrees of freedom, but additional freedom is also detrimental to the confidence in the resulting estimates due to the presence of redundancies among the parameters. Ideally, one would like to consider all possible parameter combinations and run the estimation for each one of them, yet this is computationally intractable for more than a handful of parameters. It should also be noted that an acceptable compromise may not always be found in practice. A parameter identification procedure has been developed in this work for the systematic construction of minimal parameter combinations. The idea is to start with a single parameter, and keep on adding parameters to the most promising combinations, one at a time, until some threshold on estimated parameter error is reached. To make the procedure computationally tractable, a limited number of parameter combinations are retained after each stage, for example nmax comb = 3. This selection is based on two criteria:
A formal statement of the parameter identification procedure is given in Algorithm 1 below. INITIALIZATION: • Select parameter subset, {pk , k = 1, . . . , np } rel • Select relative confidence threshold, δmax • Select max. number of parameter combinations, nmax comb • Set Pi ← ∅, i = 1, . . . , nmax comb REPEAT: • Set S ← ∅ • LOOPS: i = 1, . . . , nmax comb , k = 1, . . . , np · IF: pk ∈ Pi OR Pi ∪ {pk } has already been considered; THEN: Loop · Estimate the parameters in Pi ∪ {pk } by solving (7) · IF: Confidence criterion (8) is satisfied; THEN: Insert Pi ∪ {pk } in S • IF: S is empty; STOP • Rank parameter combination in S according to the leastsquare error criterion J in (7) • Update Pi , i = 1, . . . , min{nmax comb , |S|} with the best parameter combinations in S Algorithm 1: Systematic parameter identification procedure. 4.2 Application to Anaerobic Digestion of Microalgae The proposed parametric identification procedure is applied to the dynamic model (4)-(6) based on the experimental data set. Default values for all the parameters and initial conditions in this model are listed in Tab. 2. • Values for the stoichiometric parameters αi have been deduced from those in the ADM1 model (Batstone et al., 2002; Mairet et al., 2011); • Values for the kinetic parameters have been taken/deduced from the ADM1 model, except KI3 that is from (Bernard et al., 2001); • Initial conditions for S3 , N , Z and C have been deduced from the experimental data; regarding COD fractionation, an arbitrary value of 1 gCOD /L has been selected for the state variables S1 , S2 , X1 , X2 and X3 , based on the total COD value of 6.48 gCOD /L. • Finally, component fractions in the algal feed have been set to 30% Carbohydrates/Lipids, 40% Protein, and 30% inert COD (on a gCOD basis), based on the work of Mairet et al. (2011).
(1) Is the relative confidence of each parameter in a combination above a user-specified threshold? (2) How good is the fit in terms of the least-square error criterion? The latter is used to rank the parameter combinations and can be directly obtained from the solution of the optimization problem (7). On the other hand, the former is used to detect and eliminate those parameter combinations which yield low confidence estimates. A criterion similar to the C-optimality criteria discussed by Walter (1990) is used as a means of estimating parameter confidence. Specifically, the confidence intervals for each parameter are estimated using Student’s t-test, p Ci,i rel δi = t1−α,ν , i = 1, . . . , np , pi 2 −1 is (an estimate of) the covariance where C = Jν ddpJ2 matrix of the parameter estimate errors; ν = N × ny − np , the number of degrees of freedom; and (1 − α), the desired confidence level, e.g., 1 − α = 95%. Then, the parameter confidence criterion simply reads: rel max{δirel : i = 1, . . . , np } ≤ δmax , (8) rel with δmax the user-specified, relative confidence threshold; e.g., rel δmax = 50%.
The results reported below are obtained by considering the measurements of total gas flow rate, methane fraction, total COD, inorganic nitrogen and VFA concentrations; the pH and inert charge measurements, on the other hand, are used for validation purpose only. In applying the systematic identification procedure, the candidate parameters are µ ¯1 , KS1 , µ ¯2 , KS2 , µ ¯3 , KS3 , KI3 , KINH3 , and kco2 , along with the stoichiometric parameters. Given the lack of abundant COD information at initial time it was also deemed desirable to optimize substrate (S#,0 ) and biomass (X#,0 ) uniformly. Executing the procedure with nmax comb = 3 and relative error thresholds of 25% and 50% resulted in the optimal parameter combination presented in Tab. 3. The corresponding simulation results are plotted in Fig. 3. The good agreement between the model predictions and the experimental data, even for the data not used during the calibration, confirms that a three-reaction system is sufficient to explain the majority of the variability
5055
18th IFAC World Congress (IFAC'11) Milano (Italy) August 28 - September 2, 2011
25% confidence threshold
experimental data
8 6 4 2
0.7
0.05
0.6
VFA [gCOD /L]
10
0.8
0.04
0.03
0.02
0.5 0.4 0.3 0.2
0.01
0.1
0
0
0
1.2
0.05
8
inert charges [mol/L]
1
0.8
0.6
0.4
0.2
0
0.04 7.5 0.03
pH [-]
total COD [gCOD /L]
12
biogas flow [L/day]
50% confidence threshold
0.06
inorganic nitrogen [mol/L]
14
0.02
6.5 0.01
0 0
20
40
60
80
100
120
140
160
7
6 0
20
40
time [day]
60
80
100
time [day]
120
140
160
0
20
40
60
80
100
120
140
160
time [day]
Fig. 3. Calibration results for the 25% and 50% confidence thresholds 5. OPTIMIZATION OF ANAEROBIC DIGESTER
Table 2. Default values: stoichiometric parameters, kinetic parameters, and initial conditions α1 α2 α3 α4 α5 α6 α7 α8 α9 α10 α11 α12 S1,0 S2,0 S3,0 N0 Z0
12.5 0.006 11.5 0.04 12.5 11.5 0.083 0.04 13.1 0.006 0.19 0.12 1.0 1.0 0.011 0.022 0.038
gCOD /gCOD mol/gCOD gCOD /gCOD mol/gCOD gCOD /gCOD gCOD /gCOD mol/gCOD mol/gCOD gCOD /gCOD mol/gCOD mol/gCOD mol/gCOD gCOD /L gCOD /L gCOD /L mol/L mol/L
µ ¯1 KS1 µ ¯2 KS2 µ ¯3 KS3 KI3 KINH 3 kCO2 β1 β2 βi X1,0 X2,0 X3,0 C0 Si,0
0.175 0.45 0.2 0.3 0.2 0.15 16.384 0.0018 5.0 0.3 0.4 0.3 1.0 1.0 1.0 0.024 1.48
/day gCOD /L /day gCOD /L /day gCOD /L gCOD /L mol/L /day gCOD /L gCOD /L gCOD /L mol/L gCOD /L
in the data set (see §3.1). It is also interesting to note that increasing the parameter confidence threshold results in only a marginal improvement in model prediction capacity. The most noticeable improvements are a better fit of the tail of the inorganic nitrogen measurements and a reduction of the underestimation of the pH measurements. Table 3. Optimized parameter combinations 25% confidence threshold µ ¯1 0.369 ± 0.017 /day µ ¯2 0.0364 ± 0.0027 /day KS3 0.0711 ± 0.0145 gCOD/L KI3 1.98 ± 0.49 gCOD/L
50% confidence threshold µ ¯1 0.378 ± 0.019 /day KS2 6.43 ± 0.93 gCOD/L KS3 0.0545 ± 0.0135 gCOD/L KI3 1.53 ± 0.25 gCOD/L X#,0 0.722 ± 0.050 gCOD/L
The primary application of the calibrated grey-box model is for monitoring, control and optimization of the anaerobic digestion of microalgae. In this section, the problem of optimizing the rate of biogas (methane) production, ρCH4 , from the digester is considered, with emphasis on steady-state operation. Formally, the problem can be stated as a standard NLP with equality constraints defined by the equation system (4)-(6) at steady state: maximize: ρCH4 subject to: Steady-state model (4)-(6). An important difference with the stirred reactor used to conduct the experiments (see §2), is that the bacterial biomass is now supposed to be attached on a support for the most part, while only a small fraction is either not attached or detached by the liquid flow; consequently, only the latter is subject to the dilution effect. To account for such process heterogeneity in a simple way, a parameter, 0 ≤ θ ≤ 1, has been introduced that represents the non-attached biomass fraction (Bernard et al., 2001); that is, (θµi (·) − D)Xi = 0, i = 1, . . . , 3, with θ = 5%. Since methane has a very low solubility, it is also considered that all of the produced methane is immediately released to the atmosphere, ρCH4 = α11 µ3 (S3 , NH3 ) X3 . The decision variables in this optimization can be both the dilution rate and the feed composition. In the following case study, the dilution rate D is varied in a broad range, while the feed concentrations CODin , VFAin , Nin , Z in , and Cin are kept constant, equal to the experimental values given in §2. The model parameter values correspond to the best combination reported in Tab. 3 for a 50% confidence threshold.
5056
0.72
0.3
0.7 0.68
0.25
0.66 0.2 0.64 0.15
0.62
0.1
0.6 0.58
ρCH4 xCH4
0.05 0 0
0.5
1
1.5
0.56
30
dissolved NH3 [mol/L]
0.35
0.08
X1 X2 X3
25 20 15 10 5 0
2
2.5
3
3.5
0.3
VFA NH3
0.07
0.25
0.06 0.2 0.05 0.04
0.15
0.03 0.1 0.02 0.05 0.01 0
0
dilution rate [/day]
0.5
1
1.5
2
2.5
dilution rate [/day]
3
3.5
VFA[gCOD /L]
35 0.74
bacterial biomass [gCOD /L]
0.4
CH4 fraction [−]
CH4 production rate [mol/L/day]
18th IFAC World Congress (IFAC'11) Milano (Italy) August 28 - September 2, 2011
0 0
0.5
1
1.5
2
2.5
3
3.5
dilution rate [/day]
Fig. 4. Steady-state optimization results for various dilution rates. Left plot: methane production; Center plot: bacterial biomass; Right plot: inhibitions The results of this case study are presented in Fig. 5. The left plot shows that methane production is maximized at a dilution rate of about 2.05 day−1 (left axis). This tradeoff can be explained by the fact that a sufficiently large dilution rate is needed to bring enough COD in the digester, but too high a COD load then gives rise to various inhibitions (see below) as well as a biomass decrease due to dilution effects. Observe also that the methane fraction (xCH4 ) in the biogas remains fairly constant (left plot, right axis), giving a CH4 : CO2 ratio of about 2:1. A second auxiliary maximum for the methane production rate is observed at a higher dilution rate of about 2.95 day−1 . It is interesting to note that the maximum methane production is obtained while all three biomass species are present. The secondary maximum occurs after X2 is washed out due to too large a dilution rate (center plot, green line), meaning that all the produced methane at this operating point comes from the degradation of sugars and lipids.Then the methane flow rate quickly plummets down to zero as D keeps increasing, since all three biomass populations are eventually washed out (center plot). Based on these results, it should be clear that the model-based optimum is particularly sensitive to the biomass population growth rates. Special care must therefore be taken to determine accurate values for these parameters. Regarding inhibitions, it is found that the dissolved NH3 concentration is quite high (right plot, left axis), even at the optimum operating point. In particular, these concentrations greatly exceed the inhibition level KINH3 = 0.0018 mol/L. Therefore, regulating/optimizing the pH and alkalinity in the feed is likely to have a large influence on the methane production. The concentration of VFAs (right plot, right axis), on the other hand, remains low in comparison to the inhibition constant KI3 = 1.53 gCOD /L. 6. CONCLUSIONS AND FUTURE WORK In this paper, a two-step identification methodology has been presented and applied to the anaerobic digestion of microalgae (Chlorella vulgaris). In the first step—structural identification, the number of reactions needed to account for the variability in a given experimental data set is determined via PCA; in the second step—parametric identification, a new systematic procedure is proposed for the selection and estimation of optimal parameter combinations. A three-reaction model of the anaerobic digestion has been developed on application of this
methodology, which predicts the experimental data very well, while estimating the least possible number of parameters. Next, the calibrated model has been used to determine the steady-state operation that maximizes methane production, thereby gaining some insight about biogas production from the anaerobic digestion of microalgae. The use of multiple data sets will be considered for model identification and validation as part of future work. Moreover, the optimization of the digester operation will be investigated under transient conditions. REFERENCES Batstone, D., Keller, J., Angelidaki, R., Kalyuzhnyi, S., Pavlostathis, S., Rozzi, A., Sanders, W., Siegrist, H., and Vavilin, V. (2002). Anaerobic Digestion Model No.1 (ADM1). IWA Publishing, London, UK. Bernard, O. and Bastin, G. (2005). On the estimation of the pseudostoichiometric matrix for macroscopic mass balance modelling of biotechnological processes. Math Biosci, 193, 51–77. Bernard, O., Chachuat, B., H´elias, A., and Rodriguez, J. (2006). Can we assess the model complexity for a bioprocess? Theory and example of the anaerobic digestion process. Water Sci Technol, 53, 85–92. Bernard, O., Hadj-Sadok, Z., Dochain, Z., Genovesi, A., and Steyer, J.P. (2001). Dynamical model development and parameter identification for an anaerobic wastewater treatment process. Biotechnol Bioeng, 75(4), 424– 438. Bernard, O. and Queinnec, I. (2008). Dynamic models of biochemical processes: Properties of models. In D. Dochain (ed.), Automatic Control of Bioprocesses. Wiley-ISTE. Chisti, Y. (2008). Biodiesel from microalgae beats bioethanol. Trends Biotechnol, 3, 126–131. Dochain, D. and Vanrolleghem, P. (2001). Dynamical modelling and estimation in wastewater treatment processes. IWA Publishing, London, UK. Koster, I. and Lettinga, G. (1984). The influence of ammonium-nitrogen on the specific activity of polletized methanogenic sludge. Agricultural Wastes, 9(3), 205–216. Mairet, F., Bernard, O., Ras, M., Lardon, L., and Steyer, J.P. (2011). A dynamic model for anaerobic digestion of microalgae. In Proc 18th IFAC World Congress. Mussgnug, J., Klassen, V., Schl¨uter, A., and Kruse, O. (2011). Microalgae as substrates for fermentative biogas production in a combined biorefinery concept. J Biotechnol, 150, 51–56. Ras, M., Lardon, L., Sialve, B., Bernet, N., and Steyer, J.P. (2011). Experiental study on a coupled process of production and anaerobic digestion of chlorella vulgaris. Bioresource Technol, 102, 200–206. Sialve, B., Bernet, N., and Bernard, O. (2009). Anaerobic digestion of microalgae as a necessary step to make microalgal biodiesel sustainable. Biotechnol Adv, 27, 409–416. Walter, E.and Pronzato, L. (1990). Qualitative and quantitative experiment design for phenomenological models – A survey. Automatica, 26, 195–213. Wijffels, R. and Barbosa, M. (2010). An outlook on microalgal biofuels. Science, 329, 796–799.
5057
Anaerobic Digestion of Microalgae: Identification for Optimization and Control ⋆ Elliot Cameron ∗ Francis Mairet ∗∗ Olivier Bernard ∗∗ Monique Ras ∗∗∗ Laurent Lardon ∗∗∗ Jean-Philippe Steyer ∗∗∗ Benoˆıt Chachuat ∗,∗∗∗∗ ∗
Department of Chemical Engineering, McMaster University, Canada (Email: [email protected]) ∗∗ COMORE, INRIA Sophia-Antipolis, France (Email:{francis.mairet,olivier.bernard}@inria.fr) ∗∗∗ LBE, INRA, Narbonne, France (Email:{rasm,lardonl,jean-philippe.steyer}@supagro.inra.fr) ∗∗∗∗ Department of Chemical Engineering, Imperial College London, UK (Email: [email protected]) Abstract: Coupling an anaerobic digester to a microalgal culture is currently considered one of the most promising avenues towards the production of renewable bioenergy, either in the form of biodiesel or biogas. Accurate mathematical models are crucial tools to assess the potential of such coupled biotechnological processes and help optimize their design, operation and control. This paper focuses on the compartment of anaerobic digestion of microalgae. Using experimental data for the anaerobic digestion of Chlorella vulgaris, a grey-box model is developed that allows good prediction capabilities and retains low complexity. The proposed methodology proceeds in two steps, namely a structural and a parametric identification steps. The fitted model is then used to conduct preliminary optimization for the production of biogas from Chlorella vulgaris. The results provide some insight into the potential for bioenergy production from the digestion of microalgae and, more generally, the coupled process. Keywords: Bioprocess, anaerobic digestion, microalgae, model identification, optimization 1. INTRODUCTION Microalgae are currently considered one of the most promising feedstocks for biofuels (Chisti, 2008; Wijffels and Barbosa, 2010). But on the path to making large-scale microalgae culture sustainable, one needs to consider the management of large quantities of residual biomass along with the supply of large amounts of fertilizers. Not only does anaerobic digestion appear to be in an ideal position for addressing those challenges, but it also presents very favorable economic and energetic performance (Sialve et al., 2009; Mussgnug et al., 2011). The anaerobic digestion ecosystem is notoriously complex, involving hundreds of bacterial species. While many models of the anaerobic digestion processes have been developed since the 1970s, including the well-accepted ADM1 model (Batstone et al., 2002), there has been little direct study on the modeling of the anaerobic digestion of algal biomass. In this context, special attention must be paid to the low biodegradability of common microalgae and the large nitrogen content of algal biomass which, when converted into ammonia, may inhibit bacterial activity (Koster and Lettinga, 1984). The focus in this paper is on modeling of the anaerobic digestion of microalgae. Our objective is to develop and identify grey-box models that allow good prediction capabilities, while keeping their complexity as low as possible for use in optimization and control. ⋆ Elliot Cameron is grateful for graduate scholarship from NSERC and OGS. Support from the McMaster Advanced Control Consortium (MACC) is also gratefully acknowledged.
978-3-902661-93-7/11/$20.00 © 2011 IFAC
Modeling biotechnological processes is a tricky task, for which systematic methods and tools are still lacking. In practice, greybox models are often considered, which consist of two parts (Bernard and Queinnec, 2008): a first-principles part derived from mass-balance considerations; and a phenomenological part describing the biological reactions. The model development and identification methodology presented in this paper considers such a model structure and proceeds in two main steps. First, the minimum number of reactions needed to match the variability in a given experimental data set is determined. This allows selection of a consistent reaction scheme, which in turn imposes the key state variables. Next, the mass-balance equations and phenomenological laws are formulated, and the unknown stoichiometric/kinetic parameter values are estimated in order for the model predictions to match the experimental data. In particular, a systematic parametric identification procedure is proposed, whereby minimal parameter combinations are determined on account of both the fit quality and the estimated parameter confidence. This methodology is applied to a set of data obtained over 140 days on an experimental digester fed with Chlorella vulgaris (Ras et al., 2011), and the reaction scheme is based on the recent work by Mairet et al. (2011). An application of the calibrated model is for optimizing the production of biogas from Chlorella vulgaris. One such preliminary optimization study is conducted in this paper, thereby gaining some insight into the potential for bioenergy production from the digestion of microalgae and, by extension, from integrated biotechnological processes that couple a microalgal culture with an anaerobic digester.
5052
10.3182/20110828-6-IT-1002.03196
18th IFAC World Congress (IFAC'11) Milano (Italy) August 28 - September 2, 2011
where ξ, ξ in ∈ IRnξ are the state variables in the reactor liquid phase and in the feed stream; D, the dilution rate; K ∈ IRnξ ×nr , the stoichiometric matrix; r ∈ IRnr , the reaction rate vector; and q ∈ IRnξ , the gaseous rate vector.
−1
Dilution rate (d )
0.4 0.3 0.2 0.1
−1
Substrate additions (gCOD.L )
0
3.1 Determining Model Complexity 0
20
40
60
80
100
120
140
160
The first step consists in identifying the reaction network required to explain the variability of a set of experimental data. Determining the minimum number of reactions required to explain the set of experimental data can be performed according to the principal component analysis (PCA)-based methodology outlined in (Bernard and Bastin, 2005; Bernard et al., 2006).
6
4
2
0
0
20
40
60
80 time (d)
100
120
140
160
Fig. 1. Feed conditions of the anaerobic digester The remainder of the paper is organized as follows. The experimental setup and data are briefly described in §2. Then, the structural and parametric identification steps are detailed in §3 and §4, respectively, and the results obtained for the experimental data are presented. Finally, the fitted model is used to conduct preliminary optimization for the production of biogas from microalgae in §5. 2. EXPERIMENTAL METHODOLOGY Experimental Setup An experimental study of the anaerobic digestion of Chlorella vulgaris was performed over 150 days. This study involved the operation of a 1 L continuous-stirredtank anaerobic digester, with 0.1 L headspace, maintained at 35 ◦ C, and with no pH control. Feed Stream Characterization A detailed characterization of the feed stream can be found in (Ras et al., 2011; Mairet et al., 2011). Feed was introduced as slugs, and equal volumes of reactor medium were removed to maintain a constant liquid volume. A daily average of the dilution rate along with the substrate additions are shown in Fig. 1. The average concentrations of other relevant inlet components are provided in Tab. 1. Table 1. Feed Stream Composition Volatile fatty acids Inorganic nitrogen Inorganic Carbon Inert Cations - Anions pH
0.0 gCOD/L 0.011 mol/L 0.017 mol/L 0.017 mol/L 9.4
Measurement Techniques Readings of the following quantities were taken on an approximately daily basis: total COD (by colorimetric method); ion concentrations (by ion chromatography); biogas volume (by water displacement); biogas composition (by gas chromatography); and pH (by colorimetric method). Random samples were also selected on a less frequent basis for volatile fatty acid (VFA) determination. See (Ras et al., 2011) for more details. 3. STRUCTURAL MODEL IDENTIFICATION A general mass-balance model of the following form is considered for the system: ˙ ξ(t) = Kr (·) + D(t) ξ in (t) − ξ(t) − q (·) , (1)
The number of state variables that can be directly measured is quite limited in most biotechnological processes, and the available measurements often correspond to combinations of the state variables, such as total COD. Sometimes, measurements can also be related to the reaction rates, for example gaseous outflows of low soluble species such as methane. The analysis throughout this subsection is thus conducted for an output vector of the form: ′ C′ ξ y , (2) = y= C′′ r (·) y′′ where ny = n′y + n′′y is the total number of measurements, such ′ ′′ that ny ≥ nr . The matrices C′ ∈ Rny ×nξ and C′′ ∈ Rny ×nr represent the combinations for the states and reaction rates. For any time instants ti < tj , let u(·, ·), v(·, ·) and w(·, ·) be defined as: Z tj ′ in ′ ′ ′ D(t)[y (t) − y (t)]dt y (t ) − y (t ) − j i ti u(ti , tj ) = Z tj ′′ y (t)dt t Z tj Z tj i q(·)dt. r(·)dt, and w(ti , tj ) = v(ti , tj ) = ti
ti
Suppose that N records of u(·, ·), v(·, ·) and w(·, ·) with N > nr are available, and define the matrices U = (u(t0 , t1 ), . . . , u(tN −1 , tN )) and V, W alike. Then, from (1) and (2), the following linear relationship holds: ′ ′ CK C U= V+ W. (3) 0 C′′ | {z } =: Γ
Under the additional assumptions that the matrix Γ is full rank—which requires that the measurements y be nonredundant—and that the measured liquid species y′ are not involved in liquid-gas transfer— which enforces W = 0—the number of reactions required to describe the data can be directly be assessed from the PCA of U. Specifically, each principal axis is representative of a given reaction, and the corresponding principal component represents the relative proportion of the overall variability in the data set that can be accounted for by that reaction. Dividing by the largest eigenvalue makes it easy to determine cumulative variability, and thus the minimum number of reactions required to account for a variability in the data greater than a given threshold; e.g., 95% of the variability. In making this PCA analysis, care should always be taken to normalize and center the data in each row of U as:
5053
18th IFAC World Congress (IFAC'11) Milano (Italy) August 28 - September 2, 2011
µ1 (S1 ) = µ ¯1
100
S1 , S1 + KS1
cumulated variability [%]
90
µ3 (S3 , NH3 ) = µ ¯3
80 70
50 40 30 20 10
0
1
2
3
4
number of reactions
Fig. 2. PCA results e i = U√i − µ(Ui ) , U N σ(Ui )
i = 1, . . . , ny .
Application to Anaerobic Digestion of Microalgae The matrix U is calculated using the available liquid-phase measurements in total COD, inorganic nitrogen and VFA concentrations, along with the methane flow rate measurements that are representative of the methanogenic bacteria activity; moreover, the time instants t1 , t2 , . . . , tN +1 are directly taken as the measurement times, and cubic spline interpolants are used to compute the integral terms. Note that none of the measured liquid-phase species are involved in liquid-gas transfer, and therefore the PCA results can be interpreted in terms of biological reactions. PCA is applied to the centered and normalized data, and the cumulative variability for the principal components is presented in Fig. 2. These results indicate that a minimum of 2 or 3 reactions are sufficient for explaining, respectively, 95% and 99% of the variability in the data. For the purposes of this work a three reaction scheme is considered subsequently. 3.2 Three-Reaction Model Having settled on the number of reactions required to describe the available data, the next step is to determine what three reactions adequately perform this function. It should be noted that much individual expertise associated with grey box modeling may come into play at this stage, and this choice can be somewhat trial-and-error. In this work, the three-reaction model recently developed by Mairet et al. (2011) is considered. This model involves two hydrolysis+acetogenesis steps in parallel, whereby sugars+lipids (S1 ) and protein (S2 ) are converted into VFA (S3 ), followed by a methanogenic step. It is assumed that each reaction is performed by specific bacterial populations X1 , X2 and X3 . • Hydrolysis+acetogenesis of sugars+lipids: µ1 (·) X1
α1 S1 + α2 NH+ −→ X1 + α3 S3 + α4 CO2 4 • Hydrolysis+acetogenesis of proteins: µ2 (·) X2
α5 S2 −→ X2 + α6 S3 + • Methanogenesis: µ3 (·) X3
α7 NH+ 4
+ α8 CO2
α9 S3 + α10 NH+ −→ X3 + α11 CH+ 4 4 + α12 CO2
The biological kinetics are taken as:
S3 S3 + KS3 +
S32 KI3
(4)
Observe that knowledge of the pH is required to compute the dissolved NH3 concentration in µ3 . The following chargebalance equation assumes that all acid/base pairs are in equilibrium with h := 10−pH , and that the pH range of operation is lower than 8: h KC Z +h+ N− C h + KN KC + h KVFA− K H2 O − S3 − = 0, (5) γVFA (KVFA− + h) h with C and N denoting total inorganic carbon and nitrogen, respectively; KC , KN , KVFA and KH2 O , the dissociation con+ − stants for HCO− 3 /CO2 , NH3 /NH4 , VFA /HVFA and water, respectively; Z, the difference between inert cation and anion charges; and γVFA = 64 gCOD /mol, by assimilating VFA to pure acetate.
60
0
S2 , S2 + KS2 KINH3 . KINH3 + NH3
µ2 (S2 ) = µ ¯2
Under perfect mixing conditions, the mass-balance equations for the state variables S1 , X1 , S2 , X2 , S3 , X3 , N , C, and Z, along with the concentration Si of inert COD, read: S˙ 1 = D(β1 S in − S1 ) − α1 µ1 X1 (6) X˙ 1 = (µ1 − D)X1 S˙ 2 = D(β2 S in − S2 ) − α5 µ2 X2 X˙ 2 = (µ2 − D)X2
S˙ 3 = − D S3 + α3 µ1 X1 + α6 µ2 X2 − α9 µ3 X3 X˙ 3 = (µ3 − D)X3 N˙ = D(N in − N ) − α2 µ1 X1 + α7 µ2 X2 − α10 µ3 X3
C˙ = D(C in − C) + α4 µ1 X1 + α8 µ2 X2 + α12 µ3 X3 − ρCO2 Z˙ = D(Z in − Z) S˙ i = D(βi S in − Si ).
This formulation neglects the bacterial decay rates. It also assumes negligible NH3 and VFA gaseous losses. Expressions of the gas outflow rate and composition are given in (Mairet et al., 2011). 4. PARAMETRIC MODEL IDENTIFICATION
Once a candidate model structure has been identified, adequate values must be given to the parameters in the corresponding equation system. While accurate values for some of these parameters can be found in the literature, such as acid/base dissociation constants, other parameters may be more uncertain and need to be estimated from experimental data; this typically includes a large number of stoichiometric and kinetic parameters as well as initial conditions for the state variables. 4.1 Systematic Parameter Identification Procedure The heart of parametric model identification is to determine estimates for the uncertain parameters, which are denoted by p subsequently. Such estimates can be obtained by minimizing a specified objective function, typically a weighted sum of the square errors between measured outputs ym and corresponding model predictions y at a number of observation times tj , j = 0, . . . , N :
5054
18th IFAC World Congress (IFAC'11) Milano (Italy) August 28 - September 2, 2011
min J := p
ny N X X i=1 j=0
2 ωi,j yi (tj ; p) − yjm (tk ) ,
(7)
where ωi,j denotes the weighting coefficient for the jth measurement of the ith output; for example, these weights can be chosen as the inverse of the standard deviation of the corresponding output measurement. Moreover, the model outputs y correspond to some state variables ξ or, more generally, functions of the state variables (see §3.1). Mathematical models for biotechnological processes often comprise many uncertain parameters and relatively few measured outputs, which makes them notoriously hard to calibrate due to structural/practical identifiability issues (Dochain and Vanrolleghem, 2001). For this reason, attempting to fit all the parameters simultaneously almost always results in nearly zero confidence in the estimated parameters. Instead, two typically conflicting objectives must be accommodated in calibrating biotechnological process models, namely the fit quality and the estimated parameter confidence. Obtaining a close fit between the model predictions and the experimental data is always favored by increasing the number of degrees of freedom, but additional freedom is also detrimental to the confidence in the resulting estimates due to the presence of redundancies among the parameters. Ideally, one would like to consider all possible parameter combinations and run the estimation for each one of them, yet this is computationally intractable for more than a handful of parameters. It should also be noted that an acceptable compromise may not always be found in practice. A parameter identification procedure has been developed in this work for the systematic construction of minimal parameter combinations. The idea is to start with a single parameter, and keep on adding parameters to the most promising combinations, one at a time, until some threshold on estimated parameter error is reached. To make the procedure computationally tractable, a limited number of parameter combinations are retained after each stage, for example nmax comb = 3. This selection is based on two criteria:
A formal statement of the parameter identification procedure is given in Algorithm 1 below. INITIALIZATION: • Select parameter subset, {pk , k = 1, . . . , np } rel • Select relative confidence threshold, δmax • Select max. number of parameter combinations, nmax comb • Set Pi ← ∅, i = 1, . . . , nmax comb REPEAT: • Set S ← ∅ • LOOPS: i = 1, . . . , nmax comb , k = 1, . . . , np · IF: pk ∈ Pi OR Pi ∪ {pk } has already been considered; THEN: Loop · Estimate the parameters in Pi ∪ {pk } by solving (7) · IF: Confidence criterion (8) is satisfied; THEN: Insert Pi ∪ {pk } in S • IF: S is empty; STOP • Rank parameter combination in S according to the leastsquare error criterion J in (7) • Update Pi , i = 1, . . . , min{nmax comb , |S|} with the best parameter combinations in S Algorithm 1: Systematic parameter identification procedure. 4.2 Application to Anaerobic Digestion of Microalgae The proposed parametric identification procedure is applied to the dynamic model (4)-(6) based on the experimental data set. Default values for all the parameters and initial conditions in this model are listed in Tab. 2. • Values for the stoichiometric parameters αi have been deduced from those in the ADM1 model (Batstone et al., 2002; Mairet et al., 2011); • Values for the kinetic parameters have been taken/deduced from the ADM1 model, except KI3 that is from (Bernard et al., 2001); • Initial conditions for S3 , N , Z and C have been deduced from the experimental data; regarding COD fractionation, an arbitrary value of 1 gCOD /L has been selected for the state variables S1 , S2 , X1 , X2 and X3 , based on the total COD value of 6.48 gCOD /L. • Finally, component fractions in the algal feed have been set to 30% Carbohydrates/Lipids, 40% Protein, and 30% inert COD (on a gCOD basis), based on the work of Mairet et al. (2011).
(1) Is the relative confidence of each parameter in a combination above a user-specified threshold? (2) How good is the fit in terms of the least-square error criterion? The latter is used to rank the parameter combinations and can be directly obtained from the solution of the optimization problem (7). On the other hand, the former is used to detect and eliminate those parameter combinations which yield low confidence estimates. A criterion similar to the C-optimality criteria discussed by Walter (1990) is used as a means of estimating parameter confidence. Specifically, the confidence intervals for each parameter are estimated using Student’s t-test, p Ci,i rel δi = t1−α,ν , i = 1, . . . , np , pi 2 −1 is (an estimate of) the covariance where C = Jν ddpJ2 matrix of the parameter estimate errors; ν = N × ny − np , the number of degrees of freedom; and (1 − α), the desired confidence level, e.g., 1 − α = 95%. Then, the parameter confidence criterion simply reads: rel max{δirel : i = 1, . . . , np } ≤ δmax , (8) rel with δmax the user-specified, relative confidence threshold; e.g., rel δmax = 50%.
The results reported below are obtained by considering the measurements of total gas flow rate, methane fraction, total COD, inorganic nitrogen and VFA concentrations; the pH and inert charge measurements, on the other hand, are used for validation purpose only. In applying the systematic identification procedure, the candidate parameters are µ ¯1 , KS1 , µ ¯2 , KS2 , µ ¯3 , KS3 , KI3 , KINH3 , and kco2 , along with the stoichiometric parameters. Given the lack of abundant COD information at initial time it was also deemed desirable to optimize substrate (S#,0 ) and biomass (X#,0 ) uniformly. Executing the procedure with nmax comb = 3 and relative error thresholds of 25% and 50% resulted in the optimal parameter combination presented in Tab. 3. The corresponding simulation results are plotted in Fig. 3. The good agreement between the model predictions and the experimental data, even for the data not used during the calibration, confirms that a three-reaction system is sufficient to explain the majority of the variability
5055
18th IFAC World Congress (IFAC'11) Milano (Italy) August 28 - September 2, 2011
25% confidence threshold
experimental data
8 6 4 2
0.7
0.05
0.6
VFA [gCOD /L]
10
0.8
0.04
0.03
0.02
0.5 0.4 0.3 0.2
0.01
0.1
0
0
0
1.2
0.05
8
inert charges [mol/L]
1
0.8
0.6
0.4
0.2
0
0.04 7.5 0.03
pH [-]
total COD [gCOD /L]
12
biogas flow [L/day]
50% confidence threshold
0.06
inorganic nitrogen [mol/L]
14
0.02
6.5 0.01
0 0
20
40
60
80
100
120
140
160
7
6 0
20
40
time [day]
60
80
100
time [day]
120
140
160
0
20
40
60
80
100
120
140
160
time [day]
Fig. 3. Calibration results for the 25% and 50% confidence thresholds 5. OPTIMIZATION OF ANAEROBIC DIGESTER
Table 2. Default values: stoichiometric parameters, kinetic parameters, and initial conditions α1 α2 α3 α4 α5 α6 α7 α8 α9 α10 α11 α12 S1,0 S2,0 S3,0 N0 Z0
12.5 0.006 11.5 0.04 12.5 11.5 0.083 0.04 13.1 0.006 0.19 0.12 1.0 1.0 0.011 0.022 0.038
gCOD /gCOD mol/gCOD gCOD /gCOD mol/gCOD gCOD /gCOD gCOD /gCOD mol/gCOD mol/gCOD gCOD /gCOD mol/gCOD mol/gCOD mol/gCOD gCOD /L gCOD /L gCOD /L mol/L mol/L
µ ¯1 KS1 µ ¯2 KS2 µ ¯3 KS3 KI3 KINH 3 kCO2 β1 β2 βi X1,0 X2,0 X3,0 C0 Si,0
0.175 0.45 0.2 0.3 0.2 0.15 16.384 0.0018 5.0 0.3 0.4 0.3 1.0 1.0 1.0 0.024 1.48
/day gCOD /L /day gCOD /L /day gCOD /L gCOD /L mol/L /day gCOD /L gCOD /L gCOD /L mol/L gCOD /L
in the data set (see §3.1). It is also interesting to note that increasing the parameter confidence threshold results in only a marginal improvement in model prediction capacity. The most noticeable improvements are a better fit of the tail of the inorganic nitrogen measurements and a reduction of the underestimation of the pH measurements. Table 3. Optimized parameter combinations 25% confidence threshold µ ¯1 0.369 ± 0.017 /day µ ¯2 0.0364 ± 0.0027 /day KS3 0.0711 ± 0.0145 gCOD/L KI3 1.98 ± 0.49 gCOD/L
50% confidence threshold µ ¯1 0.378 ± 0.019 /day KS2 6.43 ± 0.93 gCOD/L KS3 0.0545 ± 0.0135 gCOD/L KI3 1.53 ± 0.25 gCOD/L X#,0 0.722 ± 0.050 gCOD/L
The primary application of the calibrated grey-box model is for monitoring, control and optimization of the anaerobic digestion of microalgae. In this section, the problem of optimizing the rate of biogas (methane) production, ρCH4 , from the digester is considered, with emphasis on steady-state operation. Formally, the problem can be stated as a standard NLP with equality constraints defined by the equation system (4)-(6) at steady state: maximize: ρCH4 subject to: Steady-state model (4)-(6). An important difference with the stirred reactor used to conduct the experiments (see §2), is that the bacterial biomass is now supposed to be attached on a support for the most part, while only a small fraction is either not attached or detached by the liquid flow; consequently, only the latter is subject to the dilution effect. To account for such process heterogeneity in a simple way, a parameter, 0 ≤ θ ≤ 1, has been introduced that represents the non-attached biomass fraction (Bernard et al., 2001); that is, (θµi (·) − D)Xi = 0, i = 1, . . . , 3, with θ = 5%. Since methane has a very low solubility, it is also considered that all of the produced methane is immediately released to the atmosphere, ρCH4 = α11 µ3 (S3 , NH3 ) X3 . The decision variables in this optimization can be both the dilution rate and the feed composition. In the following case study, the dilution rate D is varied in a broad range, while the feed concentrations CODin , VFAin , Nin , Z in , and Cin are kept constant, equal to the experimental values given in §2. The model parameter values correspond to the best combination reported in Tab. 3 for a 50% confidence threshold.
5056
0.72
0.3
0.7 0.68
0.25
0.66 0.2 0.64 0.15
0.62
0.1
0.6 0.58
ρCH4 xCH4
0.05 0 0
0.5
1
1.5
0.56
30
dissolved NH3 [mol/L]
0.35
0.08
X1 X2 X3
25 20 15 10 5 0
2
2.5
3
3.5
0.3
VFA NH3
0.07
0.25
0.06 0.2 0.05 0.04
0.15
0.03 0.1 0.02 0.05 0.01 0
0
dilution rate [/day]
0.5
1
1.5
2
2.5
dilution rate [/day]
3
3.5
VFA[gCOD /L]
35 0.74
bacterial biomass [gCOD /L]
0.4
CH4 fraction [−]
CH4 production rate [mol/L/day]
18th IFAC World Congress (IFAC'11) Milano (Italy) August 28 - September 2, 2011
0 0
0.5
1
1.5
2
2.5
3
3.5
dilution rate [/day]
Fig. 4. Steady-state optimization results for various dilution rates. Left plot: methane production; Center plot: bacterial biomass; Right plot: inhibitions The results of this case study are presented in Fig. 5. The left plot shows that methane production is maximized at a dilution rate of about 2.05 day−1 (left axis). This tradeoff can be explained by the fact that a sufficiently large dilution rate is needed to bring enough COD in the digester, but too high a COD load then gives rise to various inhibitions (see below) as well as a biomass decrease due to dilution effects. Observe also that the methane fraction (xCH4 ) in the biogas remains fairly constant (left plot, right axis), giving a CH4 : CO2 ratio of about 2:1. A second auxiliary maximum for the methane production rate is observed at a higher dilution rate of about 2.95 day−1 . It is interesting to note that the maximum methane production is obtained while all three biomass species are present. The secondary maximum occurs after X2 is washed out due to too large a dilution rate (center plot, green line), meaning that all the produced methane at this operating point comes from the degradation of sugars and lipids.Then the methane flow rate quickly plummets down to zero as D keeps increasing, since all three biomass populations are eventually washed out (center plot). Based on these results, it should be clear that the model-based optimum is particularly sensitive to the biomass population growth rates. Special care must therefore be taken to determine accurate values for these parameters. Regarding inhibitions, it is found that the dissolved NH3 concentration is quite high (right plot, left axis), even at the optimum operating point. In particular, these concentrations greatly exceed the inhibition level KINH3 = 0.0018 mol/L. Therefore, regulating/optimizing the pH and alkalinity in the feed is likely to have a large influence on the methane production. The concentration of VFAs (right plot, right axis), on the other hand, remains low in comparison to the inhibition constant KI3 = 1.53 gCOD /L. 6. CONCLUSIONS AND FUTURE WORK In this paper, a two-step identification methodology has been presented and applied to the anaerobic digestion of microalgae (Chlorella vulgaris). In the first step—structural identification, the number of reactions needed to account for the variability in a given experimental data set is determined via PCA; in the second step—parametric identification, a new systematic procedure is proposed for the selection and estimation of optimal parameter combinations. A three-reaction model of the anaerobic digestion has been developed on application of this
methodology, which predicts the experimental data very well, while estimating the least possible number of parameters. Next, the calibrated model has been used to determine the steady-state operation that maximizes methane production, thereby gaining some insight about biogas production from the anaerobic digestion of microalgae. The use of multiple data sets will be considered for model identification and validation as part of future work. Moreover, the optimization of the digester operation will be investigated under transient conditions. REFERENCES Batstone, D., Keller, J., Angelidaki, R., Kalyuzhnyi, S., Pavlostathis, S., Rozzi, A., Sanders, W., Siegrist, H., and Vavilin, V. (2002). Anaerobic Digestion Model No.1 (ADM1). IWA Publishing, London, UK. Bernard, O. and Bastin, G. (2005). On the estimation of the pseudostoichiometric matrix for macroscopic mass balance modelling of biotechnological processes. Math Biosci, 193, 51–77. Bernard, O., Chachuat, B., H´elias, A., and Rodriguez, J. (2006). Can we assess the model complexity for a bioprocess? Theory and example of the anaerobic digestion process. Water Sci Technol, 53, 85–92. Bernard, O., Hadj-Sadok, Z., Dochain, Z., Genovesi, A., and Steyer, J.P. (2001). Dynamical model development and parameter identification for an anaerobic wastewater treatment process. Biotechnol Bioeng, 75(4), 424– 438. Bernard, O. and Queinnec, I. (2008). Dynamic models of biochemical processes: Properties of models. In D. Dochain (ed.), Automatic Control of Bioprocesses. Wiley-ISTE. Chisti, Y. (2008). Biodiesel from microalgae beats bioethanol. Trends Biotechnol, 3, 126–131. Dochain, D. and Vanrolleghem, P. (2001). Dynamical modelling and estimation in wastewater treatment processes. IWA Publishing, London, UK. Koster, I. and Lettinga, G. (1984). The influence of ammonium-nitrogen on the specific activity of polletized methanogenic sludge. Agricultural Wastes, 9(3), 205–216. Mairet, F., Bernard, O., Ras, M., Lardon, L., and Steyer, J.P. (2011). A dynamic model for anaerobic digestion of microalgae. In Proc 18th IFAC World Congress. Mussgnug, J., Klassen, V., Schl¨uter, A., and Kruse, O. (2011). Microalgae as substrates for fermentative biogas production in a combined biorefinery concept. J Biotechnol, 150, 51–56. Ras, M., Lardon, L., Sialve, B., Bernet, N., and Steyer, J.P. (2011). Experiental study on a coupled process of production and anaerobic digestion of chlorella vulgaris. Bioresource Technol, 102, 200–206. Sialve, B., Bernet, N., and Bernard, O. (2009). Anaerobic digestion of microalgae as a necessary step to make microalgal biodiesel sustainable. Biotechnol Adv, 27, 409–416. Walter, E.and Pronzato, L. (1990). Qualitative and quantitative experiment design for phenomenological models – A survey. Automatica, 26, 195–213. Wijffels, R. and Barbosa, M. (2010). An outlook on microalgal biofuels. Science, 329, 796–799.
5057