Can the operating limits of biogas plants operated under non-isothermal conditions be defined with certainty? Modeling self-optimizing attainable regions.

ContentslistsavailableatScienceDirect

Computers and Chemical Engineering

journalhomepage:www.elsevier.com/locate/compchemeng

Can the operating limits of biogas plants operated under

non-isothermal conditions be defined with certainty? Modeling self-optimizing attainable regions

F. Abunde Neba

^a,d,e,∗

, Michel Tornyeviadzi

^a,e

, Nana Y. Asiedu

^c

, Ahmad Addo

^b

, John Morken

^f

, Stein W. Østerhus

^d

, Razak Seidu

^e

aAbunde Sustainable Engineering Group, (AbundeSEG), Cameroon

bDepartment of Agricultural and Biosystems Engineering, Kwame Nkrumah University of Science and Technology, Kumasi, Ghana

cDepartment of Chemical Engineering, Kwame Nkrumah University of Science and Technology, Kumasi, Ghana

dDepartment of Civil and Environmental Engineering, Norwegian University of Science and Technology, Trondheim, Norway

eDepartment of Marine Operations and Civil Engineering, Norwegian University of Science and Technology, ˚Alesund, Norway

fFaculty of Science and Technology, Drobakveien 31, 1432 Aas, Norwegian University of Life Sciences, ˚As, Norway

a rt i c l e i n f o

Article history:

Received 18 October 2019 Revised 31 January 2020 Accepted 29 June 2020 Available online 30 June 2020 KeyWords:

Self-optimizing operating limits Attainable regions

Process kinetics Biogas digesters Temperature

a b s t r a c t

Uncertaintyinoperatingparameterssuchastemperatureunderminesthereliabilityofusingkineticmod- elsinperformanceprojectionsforplantsoperatedunderambientnon-isothermalconditions.Thisstudy developsatheoreticalframework,whichusesprocesskinetics,uncertaintyquantificationtodefinerobust operatinglimits knownasself-optimizingattainableregions,where byinsteadofdefiningaverylarge operatinglimit,whichwillbe achievedsomeofthetimesfor someofthe reactorconfigurations,we defineaself-optimizinglimit,whichwillbeachievedallthetimesforallpossiblereactorconfigurations (despitevariationsintemperature).Usingatemperaturerangeof20– 60^◦C,,theresultsindicatethatde- creasingtemperatureuncertainty,increasingprocesstemperatureorusingamultistagedigesterstructure increasestheself-optimizing operatinglimits:1.53×10⁻⁴,4.95×10⁻⁴ and6.32×10⁻⁴ (g/L)^{2 obtained} fortemperaturesof20.00,31.60and52.40^◦Crespectively.Thefindingshighlyimportantindefiningper- formancetargetsespeciallywhenthereisuncertaintyinenvironmentalconditions.

© 2020 The Authors. Published by Elsevier Ltd.

ThisisanopenaccessarticleundertheCCBYlicense.(http://creativecommons.org/licenses/by/4.0/)

1. Introduction

Environmental pollution resulting from the use of fossil fuels coupled with potential exhaustion of fossil fuel resources makes itnecessary tofindalternativeenergysources thatare renewable andenvironmentallysustainable.Biogas,amethane-richrenewable energy obtained fromthe anaerobic treatment of organic wastes hasproventohaveagreat potentialinshiftingrelianceonfossil- basedenergy.However,biogasproductionfromanaerobicdigestion isinfluencedby amyriadoffactorsincludingsubstratecharacter- isticsandconcentration,presenceofinhibitorysubstancesaswell as temperature and pH and reactor configurations (Henze et al., 2008; Wang et al., 2007). Amongst these factors, temperatureis consideredveryimportant,sinceitinfluencestheactivityofanaer- obic microorganisms by influencing the activity of some essen-

∗ Corresponding author.

E-mail address: [email protected] (F. Abunde Neba).

tial enzymes involved in organic matter degradation and biogas production (Kim et al., 2017; Donoso-Bravo et al., 2013). Biogas plantscan generallybe operated underfourtemperatureregimes dependingon the type of microorganisms present: psychrophilic digestion(15–25^◦C), mesophilic digestion(30–40^◦C), thermophilic digestion (50–60 ^◦C) and hyper-thermophilic (65–75 °C) diges- tion (Kuo and Lai, 2010; Wang et al., 2012; Saady and Massé, 2015; Pavlostathis andGiraldo-Gomez, 1991). Within each of the temperatureregimes,theactivityofanaerobicmicroorganismsin- creases with increasing temperature up to a maximum activity abovewhichasharpdropisobservedwithincreaseintemperature (PavlostathisandGiraldo-Gomez,1991). Intermsofcapacity,bio- gasplantscanusually beoperatedaslarge-scale(common inthe developedworld)orsmall-scale(commonindevelopingcountries) systems.Thelarge-scalesystemsarenormallyoperatedundercon- trolled isothermal conditions (most commonly at mesophilic or thermophilicregimes).Thisimpliesthat theinfluenceofthetem- peratureonlarge-scaleplantscanbeconsideredconstant(Donoso- Bravoetal.,2013),andhencethe performanceoflarge-scale sys- https://doi.org/10.1016/j.compchemeng.2020.107001

0098-1354/© 2020 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license. ( http://creativecommons.org/licenses/by/4.0/ )

Nomenclature

(VS)₀ Initialconcentrationofvolatilesolids(gVS/L) A_f Acidityfactor(gVFA/L)/(gBVS/L)

B₀ Biodegradabilityconstant(gBVS/L)/(gVS/L) K_ime VFA inhibition constant for methanogenic archae

(gVFA/L)

K_sac Monod half-saturation constant for acidogenic bacteria(gBVS/L)

Ks_me Monod half-saturation constant for acidogenic bacteria(gVFA/L)

S_BVS

0 Initial concentration of biodegradable volatile solids(gBVS/L)

SBVS Concentration of biodegradable volatile solids (gBVS/L)

S_VFA

0 Initialconcentrationofvolatilefattyacids(gVFA/L) S_VFA Concentrationof volatilefattyacids inbioreactor

(gVFA/L)

Tmax Maximum temperature at which growth rate is zero(^◦C)

T_min Minimum temperature at which growth rate is zero(^◦C)

X₀ Initialconcentrationofbiomassinreactor(g/L) X_ac0 Initialconcentrationofacidogenicbacteria(gac./L) Xac Concentrationofacidogenicbacteriainbioreactor

(gac./L)

Xme₀ Initial concentration of methanogenic archaea (gme./L)

Xme Concentration ofmethanogenic archaea in biore- actor(gme./L)

k₁ Yieldconstant(gBVS/gac./L) k₂ Yieldconstant(gVFA/gac./L) k₃ Yieldconstant(gVFA/gme./L) r_S

BV S Reaction rate for biodegradable volatile solids (gBVS/L/d)

r_{SV FA} Reactionrateforvolatilefattyacids(gVFA/L/d) r_Xac Reactionrateforacidogenicbacteria(gac./L/d) r_Xme Reactionrateformethanogenicarchae(gme./L/d) t_υ,α/2 Studentt-distributionparameter

β

ˆ ^{Vectorofestimatedmodelparameters}

γ

CH₄ Volumetricmethaneproductivity(LCH₄/m³/d)

γ

^{s Methaneyield}

μ

mac Maximumspecificgrowthrateofacidogenicbac- teria(d⁻¹)

μ

^mme Maximum specific growth rate of methanogenic archaea(d⁻¹)

μ

^{ac Specific growth rate of} methanogenic archaea (d⁻¹)

μ

me Specific growth rate of methanogenic archaea (d⁻¹)

σ

^{2 Standarderror}

B Ratkowskyparameter(^{◦C−1 h−0.5}) C RatkowskyparameterC(^◦C−1) T Reactortemperature(^◦C)

EMY₉₀ 90%experimentalmethaneyield(mLCH₄/gVS) HRT(T₉₀) Hydraulic retention time that gives 90% experi-

mentalmethaneyield(d)

J Jacobianmatrixevaluatedatparameterestimates VSL Volatilesolidsloading(gVS/l)

VSR Volatilesolidsreduction(%) n Numberofexperimentaldatapoints p Numberofmodelparameters

α

Significancelevel

β

^{Vectorofrealmodelparameters}

ϑ Acidogenicfraction

tems can more closely mimictheoretical predictions duringeco- nomic feasibilitystudies,which useprojectionsinbiogasproduc- tion corresponding to a fixed temperature. On the other hand, small scale systems (such as domestic biogas septic tanks) are mostlyoperatedunderambientnon-isothermalconditions(mostly duetocostconstraints),whereheatingofhighvolumesoforganic wasteisunlikely(Bandaraetal.,2012;Suminoetal.,2007).Small scale plants operated under these non-isothermal conditions are usually viable to uncertain performance resulting from tempera- turefluctuationsmainlyduetoseasonalvariations,day–nightpat- terns as well as specific spikes or drops events. This makes the definition of operating limits a challenging task for small scale systemsaslargeperformance deviationsare usuallyobserved be- tween predictedand actual. There have been several studies in- vestigating the influence of temperature on the performance of the anaerobictreatment process examining microbial community (Kimetal.,2017),organicmatterremoval(Lohanietal.,2018),bio- gasproduction(Lietal., 2013;Latifetal., 2012)anddigestercon- figuration(Chongetal.,2012).Thestudies investigatingtheeffect oftemperature onthe performance ofanaerobicdigestion gener- allyproceed by fourmain steps:(1) Defininga digester configu- ration or listof digester configurations to be tested, (2) defining alistoftemperatureswithinthepsychrophilic,mesophilicorther- mophilicrange(3)evaluatingtheprocessperformance(e.g.organic matterremoval,biogasproductionornutrientrecovery)forthede- finedtemperaturesanddigesterconfigurations(4)Makingadeci- sionon whichtemperaturegivesthe bestperformance forwhich digesterconfiguration.Theresultsofsuchstudiesareonlyrelevant to industrial scale projections wheretemperature is usually con- trolled,andtheoptimalsolutionismostprobablyalocaloptimum since it is selected from the options predefined insteps (1)and (2).What ifother digesterconfigurations existthat couldperform better fortemperatures not considered?Severalconfigurations of anaerobicdigesters havebeen developedaimed atenhancing the performanceofthetreatmentprocess withregardstobiogaspro- duction ororganicmatter removal. Ageneralizedclassification of anaerobic digesters configurations fall under three main groups:

single-stagesystems,multi-stagesystems,andhybridsystemssuch aselectrochemicalanaerobicdigesters(DeVriezeetal.,2018),solar bio-hybriddigesters(BustamanteandLiao,2017),etc.Amongstthe aforementioned configurations, severalstudies have reported the multiple-staged systems to be the mostpromising technology in termsofflexibilityandabilitytooptimizeeverystepoftheanaer- obic digestionprocess (Akobi etal., 2016; Zhang et al., 2017). In multistage digestion systems, the concept is to stage the anaer- obic digestion process based on kinetic and physiological differ- ences betweenthe differentgroups ofmicroorganisms catalysing thedifferentstagesof theAD process(EPA,2006). Mostresearch on staged anaerobic digestion has focused on comparative anal- ysis betweensingle and multiple stage anaerobic digesters, with less effort geared towards obtaining an optimal staging configu- ration.Whiletheperformance advantagesofmultistageanaerobic systems as opposedto single staged havebeen well established, thereislackofa systematicdigester designapproachesthatcon- siders the temperaturevariations on the performance multistage systems.Currentapproachesareempirical,involvingexperimental evaluationof predefinedconfigurations of differenttemperatures, which is not only time consume and expensive, but also results in a local optimum. An approach to determine the performance targetsforall possibletemperaturevariationsandforall possible digester configurations will therefore be a breakthrough to sup-

port investment decisions in biogas plants operated under non- isothermalconditions.

Unlike conventional studies that would proceed by evaluat- ing the effect of temperature on a given digester configuration, we rather proceed by defining the performance targets for all possible digester configurations and for all possible temperature fluctuations. The current study discusses how the attainable re- gion concept could be used for performance targeting and di- gester synthesis under non-isothermal plant operation. The at- tainable region (AR) represents a collection of all possible out- puts for all possible process configurations by interpreting pro- cesses asgeometricobjectsthat define theoperatinglimitswith- out having to explicitly enumerate all possible design configu- rations (Hildebrandt and Glasser, 1990; Hildebrandt et al., 1990; Ming et al., 2013; Ming et al., 2016). Our recent studies, have beenfirst oftheir kindillustrating theusefulnessofAR tomodel digester configurations that optimize methane productivity and volatilesolidsreduction(AbundeNebaetal.,2019c)aswellassta- bilityof methanogenic archaea(AbundeNeba etal., 2019b). Both studies put together haveillustrated that a change inthe kinetic modelstructureorvalue ofkinetic coefficients,induced bydiffer- encesinsubstrateandinoculumcharacteristicssignificantly influ- encestheperformancetarget aswell astheoptimaldigestercon- figuration requiredto achievethe target.In anotherrecentstudy, the authors demonstrated the coupling ofattainable regions and fuzzy multicriteria decisionmaking forsynthesis ofhighrate di- gesterswithoutanyneedforakineticmodel(AbundeNebaetal., 2019d). Also, a framework for coupling simplified kinetic models andeconomicfeasibilityindicatorsforsynthesis ofdigesterstruc- turesusingattainableregionshavebeendemonstratedbytheau- thors (AbundeNeba et al., 2019a). The studies were able to de- fine the absolute best performance for all possible digester con- figurations considering mixing and reaction as fundamentalpro- cesses occurring in the digester. However the studies assumed isothermal conditions where the attainable region are definedat agivenoperatingtemperatureandisthereforewell suitedforin- dustrial scale operation where temperature is controlled. In the current study, we seekextend the previous works by presenting thelimitsofachievabilityfornon-isothermalbiogasplantsi.e.bio- gas plantsoperated underthermaluncertainty. Weintroduce the concept of self-optimizingoperation, an attainableregion or per- formance target that results in near optimal operation in spite of temperature fluctuationswithin the process. A self-optimizing operation (Gausemeier et al., 2006; Permin etal., 2016) is when we canachieve anacceptable lossperformanceby usingconstant setpoint valuesforoperating parameters (e.g.temperature, kinet- ics, substratecharacteristics,etc.,) withouttheneedtoreoptimize whenvariationsoccur.SincetheARtheoryinvolvesmixingandat- tainability ofstatesunderdefinedconditionsoftemperature,feed composition and kinetics,theidea behind thisstudyis tomodel the uncertainty in state predictions resulting from temperature fluctuationsandpropagatethisuncertaintyontotheattainablere- gions to define the self-optimizing attainable regions. Unlikethe attainable regions, which are valid for definedtemperatures, the self-optimizingattainableregions(thoughresultsinacceptableloss inperformance)isvalidforallpossibletemperatureswithinade- finedrage.

2. Materialsandmethods

Thisstudyappliesaproactiveapproachofdealingwiththermal uncertainty,wherebybeforedefiningthelimitsofachievabilityof thesystem(constructing theattainableregions),we firstquantify the model prediction uncertaintyresulting from thetemperature variations.Threedifferentuncertaintybounds(the10thpercentile, mean and the 90th percentile) are used to construct three sets

ofattainableregionsandtheregionofintersectionbetweenthese threesets representsthe limitsofachievability (performancetar- gets) of the system under non-isothermal conditions. Finally we apply the necessary conditions of attainable regions to interpret theboundaries ofthe intersectionregion intodigester configura- tions.Tobetter illustratethe methods,givenan organicsubstrate tobe used asfeedstockforbiogas production,Fig.1 outlines the workflowneededtodefinethemaximumdesignandperformance targetsforagivenanaerobictreatmentprocess.

The following steps provide more detailed description and/or rationalbehindeachofthestepspresentedinthemethodological workflow(Fig.1).

2.1. Organicsubstrateandcharacteristics

Thetargetsofbiogasproductiondependsonboththeplantde- signandthefeedstockcharacteristics.Thescopeofthisstudyison plantdesignandconsidersthatthesubstrateaswellasitscharac- teristics and potential to produce biogas are known. The organic substrate considered for the anaerobic digestion process is diary manure,whichwasdigestedwiththefollowingexperimentalchar- acterization (Kafle and Chen, 2016): HRT (^T90)=28days, EMY = 204mL/gVS,90% EMY=18.6mL/gVS,VSL=3.5gVS/l,VSR=58.6%. The datawas utilizedto calculated the volumetric methanepro- ductionrate(

γ

CH₄ =2.95l/m³/d)requiredforthemodelingframe- workpresentedinthisstudy.Theformulaeappliedforcalculations areaspresentedinourpreviousstudy(AbundeNebaetal.,2019c).

γ

^CH4=EMY90

HRT ×VSL (1)

2.2.Modelingtheanaerobictreatmentprocess

Thegeometricoptimizationtechniqueofattainableregionsuti- lizedin this study isbased on the process kinetics. Because the geometriccalculations(suchasratevectors,CSTRlocus,PFRtrajec- tories)arecomplex(willbeexplainedfurtherinSection3.2.3),itis requiredtouseasimplifiedmodeloftheanaerobictreatmentpro- cess,which presenta compromise betweenbeinghighlyaccurate butverycomplexinputrequirementandhighlysimplifiedbutvery limitedpredictiveability.Inaddition,becausewearedealingwith non-isothermal conditions,thenormalArrhenius equation widely usetomodel temperaturedependence willnot beappropriate as itonlyreproducestheincreasingpartofthetemperatureinfluence ontheactivityofanaerobicmicroorganismsmeanwhile theactiv- ityofmicroorganismsdropsignificantlyafterreachinganoptimum (PavlostathisandGiraldo-Gomez,1991).Therefore,selectionofap- propriatekineticandtemperaturemodelstructuresisimportant.

2.3.Modelingtemperaturesensitivity

Modeling temperature sensitivity is required to provide dy- namicinformation on howthe states ofthe anaerobictreatment process respond to temperature fluctuations. This information is usefulin identifyingtime intervalswherethe ADprocess ismost sensitivetosuchfluctuationsorwheretemperaturecontrolofthe digesteris moreor lessimportant.Forinstance,ifthe sensitivity ofmethanogenic archae to temperature is closeto zeroin some timeinterval,changesinthevalueoftemperatureatthattimein- tervalwouldhavelittleimpactonbiogasproduction.Thisinforma- tionhassignificanteconomicimportancebecauseinsteadofheat- ing the digester throughout the entire process (which is energy consuming),intermittentheatingcanbeappliedwherethesystem isonlyheatedattimeswerethestatesaremostsensitive.

Fig. 1. Methodological work flow of the study.

2.4.Uncertaintyquantificationonpredictionoutputs

The uncertaintyassociated with the anaerobic treatment pro- cessarisesfromtemperature fluctuationsmainly duetoseasonal variations, day–nightpatterns as well asspecific spikes ordrops events.Thepropagationof thisuncertaintyaswell ascalculation ofuncertaintybounds(10th,percentile,meanand90thpercentile) onthe model stateswas done using the MonteCarlo simulation method.This procedure is highlydependent on the temperature uncertainty range (known as variance metric of input).This vari- ancemetricswasobtainedbymeasuringtheambienttemperature over a two month periodand usingthe data to fit a probability distributionfromwherethetemperaturewassampledduringev- eryiterationoftheMonteCarloprocedure.

2.5.Maximumdesignandperformancetargets

Foragivenorganicsubstrateanditsconcentration,theattain- ableregionistheconvexhullofthesetofpointsachievablecon- sideringmixingandreaction (biodegradation) asthe only funda- mentalprocesses occurringintheanaerobicdigester. Thisdefines thelimitsofachievability(performancetargetsoftheprocess)and theboundaryortheARcanbeinterpretedintodigesterconfigura- tions.Itshouldbenotedthatifmixingandbiodegradationarethe onlyfundamentalprocesses occurringinthe digester,there isno needtodeviseneworperhapsnoveldigesters(otherthanthefun- damentalPFR andCSTR configurations) that can serveto expand theoperatinglimitsofthesystem(Mingetal.,2016).Byexpanding

the performance targets, the authors mean achievingstates that were initially not achievableconsidering mixingand biodegrada- tion. Therefore,we call thispoint the maximumdesign andper- formancetarget.

3. Theoreticaldevelopments

Thetheoreticaldevelopmentsisdividedintotwomainsections:

thefirstsection(Section3.1)presentsthekineticandtemperature dependencemodelsofthe anaerobictreatmentprocess whilethe section(Section3.2)presentstheuseofattainableregionanalysis todefinethemaximumdesignandperformancetargetsoftheAD process

3.1. Modelofanaerobicbiogasreactor

3.1.1. Kineticmodel

The simplified dynamic model for anaerobic digestion of an- imal manure, which has been presented in our previous study (AbundeNebaetal., 2019c) formodeling configurationsofanaer- obicdigestersusingattainableregionswasadoptedinthecurrent work.Fig.2presentsasummaryofthemodelscheme,clearlyout- liningthemodelinputs(temperature,organicload,retentiontime, waste andbiomass characteristics) model outputs (volatile solids reductionandvolumetricmethaneproductivity),kineticconstants aswell asstatevariables.Thismodelisadvantageousbecausethe geometric calculations involved in attainable region analysis are relatively complex and it is required to havea simplified model

Fig. 2. Model of methane bioreactor showing inputs, outputs, parameters and state variables ( Abunde Neba et al., 2019c ).

(without compromise for process information) in order to make theproblemmoretractable.

The model considers four state variables, which include:

Biodegradable Volatile Solids (BVS); Volatile Fatty Acids (VFA);

acidogenic bacteria (Xac) and methanogenic bacteria (Xme). The rateexpressions forthe fourstates are respectivelypresented by Eqs.(2)to(5).

dS_BVS

dt =rBVS=−k1

μ

^{acXac (2)}

dSVFA

dt =rVFA= k2

μ

acXac−k3

μ

meXme (3)

dXac

dt =rac=

μ

^{acXac (4)}

dXme

dt =rme=

μ

meXme (5)

The model assumesthe specific death rates ofboth microbial populationsarenegligiblecomparedtotheirspecificgrowthrates.

The specific growthrate of acidogenicbacteria is modeled using the Monod equation, Eq. (6) while an uncompetitive inhibition termisaddedtothatofmethanogenicbacteria,Eq.(7)toaccount forvolatileacidinhibitionduringreactorupsetorfailure.

μ

ac=

μ

mac

SBVS

Ksac+SBVS

(6)

μ

me=

μ

mme

SVFA

Ksme+S_VFA

1+^S_K^{V FA}

ime

⁽⁷⁾

The percentage volatile solids reduction and the volumetric methaneproductivity(model outputs)arecalculatedaspresented in Abunde et al. (Abunde Neba et al., 2019c). The following pa- rametersk=[k₁, k₂, k₃,

ϑ

,

γ

^{s]aretobeestimatedwhiletheval-}

uesofallotherparametersparticularlyKsac, KsmeandK_imearemain- tainedasintheoriginalHillmodel.Theparameterestimationcon- sistedofiterativelysearchingforparametervaluesthatminimizes the squarederrorbetweentheoutputs predictedby the parame- terizedmodelandobservedexperimentally,Eq.(8).

mink

ε (

^x,k

)

⁼

γ

^CH4

(

^k

)

⁻

γ

C^eH4

2

+

(

^VSR

(

^k

)

^−VSRe

)

^{2 (8)}

Forthispurpose,theMatlaboptimizationroutine,fminconwas used,wherethedynamicmethanebioreactormodelintegratednu- mericallyusingtheRunge-Kutta4–5thordermethodimplemented bytheMatlabode45routine. Table1presentstheparameteresti- matesaswellasmodelidentificationerrorcomputedforanaerobic digestionof diary manure.The calculated experimental values of model outputs, volatile solidsreduction andvolumetric methane productivity are 58.62% and 22.95 respectively while the model predictedvaluesobtainedare58.60%and22.91.

The biodegradability andacidity constantsare unique to each type of waste Hill, 1983) andtherefore serve to characterize the wastetypebeingdigested.Forthisstudy,respectivevaluesof0.40 and0.05wereusedforthebiodegradabilityandacidityconstants.

Theeffect ofbothconstantson theinitial concentration ofstates ismodeledbyEqs.(9)and((10).

SBVS0=B0

(

^VS

)

0 (9)

SVFA0=AfSBVS0 (10)

Characterizingtherawwasteby usingBo andA_fsimplifies the modelsincetheeffectsofotherparameterssuchascationconcen- tration,alkalinity, ammoniadissolved carbon dioxideare intrinsi- callypartofBoandA_f.

3.1.2. Temperaturemodel

Theeffectofthetemperatureontheanaerobictreatmentpro- cess has been modelled using the Ratkowsky expanded square root model, Eq. (11) (Ratkowsky et al., 1983). The advantage of thismodelis thatit can describe thetemperatureinfluence over the entire biokineticrange of the anaerobic digestion process as opposed to the Arrhenius model, which only reproduces the in- creasingpart of temperature dependence. The maximum growth ratesofacid-formingbacteria(

μ

^mac) andmethaneformingarchae (

μ

mme)arefunctionsofthedigestiontemperatureandthisdepen- dencewasmodeledasshowninEq.(11)

μ

mac

(

^T

)

⁼

μ

mme

(

^T

)

^=[B

(

^T−Tmin

)

^]2

{

^1−exp[C

(

^T−Tmax

)

^]

}

²

(11)

Tmin<T<Tmax

T_min and Tmax are respectively the maximum and minimum temperaturesatwhichthegrowthrateiszerowhiletheconstants

Table 1

Estimated kinetic coefficients and model identification error.

Parameter k 1 k 2 k 3 ϑ γ^{s K sac K me K} ime Model error ( ɛ ) Value 1.096 0.096 5.351 0.519 0.503 9.0 2.0 6.0 0.0435

B(^{◦C−1 h−0.5}) ^andC (^◦C−1)^{areknown asRatkowsky}parameters, whicharenormallyestimatedfromtestdatatoreflecttheprocess beingmodelled.

The determination ofthe Ratkowsky parameters (B and C) as well asT_min and T_maxwasmade inour previous study by fitting theChenandHashimoto curve,citedby Hill(Hill, 1983)fortem- peraturedependence ongrowth ratetothe Ratkowsky expanded squarerootmodel.The parametersobtainedwere0.02,0.05, 4.22 and79.96respectivelyforB,C,T_minandTmaxrespectively.

3.2.Definingmaximumdesignandperformancetargets

3.2.1. Temperaturesensitivitymodeling

Modelingtemperaturesensitivityconsistofanalyzingthesensi- tivityofthemodelstatestothetemperature.Giventhekineticand temperaturemodelfortheanaerobicdigestionprocess,thefollow- ing section illustrates how to model the temperature sensitivity ontheanaerobicdigester. Since wedo not havean explicitsolu- tiontothecompleteprocessmodel,theabsolutesensitivitiesmust becomputedusingthesensitivityequations.Forann-dimensional systemgivenbyEq.(12)

Y˙ = f

(

^t,Y;

β

, T

)

, Y

(

⁰

)

=Y₀ (12) WithY∈Rⁿ,statevariable,

β

∈R^pthemodelparameters,T∈R thetemperatureandY₀theinitialcondition,thevectoroftemper- aturesensitivities

∂

^Y/

∂

^Tsatisfy

d dt

∂

^Y

∂

^{T =}

∂

^F

∂

^Y

∂

^Y

∂

^{T +}

∂

^F

∂

^{T (13)}

Withinitialconditions

∂

^Y

(

⁰

)

∂

^{T =0n×1 (14)}

∂

^Y/

∂

^{Tisthechangeofstateswithrespectto}temperature.The sensitivityequationsarecoupledwiththeoriginalmodeldifferen- tialequationsandsolvedtoobtainthetemperaturesensitivitiesfor thenecessarytimepoints. Theresultingmatrixofabsolutesensi- tivitiesattimepointtSa(^t)=

∂

^Y/

∂

^{T willbeoftheformshownby}

Eq.(15).

Sa

(

^t

)

⁼

⎡

⎢ ⎢

⎣

S1,t1 S2,t1 S3,t1 S4,t1

S_1,t2 S_2,t2 S_3,t2 S_4,t1 ..

. ... ... ... S1,tn S2,tn S3,tn S4,tn

⎤

⎥ ⎥

⎦

⁽¹⁵⁾

3.2.2. Propagationoftemperatureuncertainty

Theglobaloptimizationtechniqueofattainableregionsisbased onageometricrepresentation(orconvexhull)ofmodelstates(or setof points) that are achievable by the system. Any factor that influences the states being output by the digester will therefore influencethe operatinglimitsofthe system, which isdefinedby theattainableregion.Thepropagationoftemperatureuncertainty onto the model stateswas done using the MonteCarlo simula- tionprocedure,whichinvolvesthreesteps:(1)specifyingtemper- ature uncertainty range, usually in the form of a statistical dis- tribution (2)sampling a definednumber ofvalues form the dis- tribution in which case we used 1000 and (3) propagating the sampledinput uncertainty on to the model states.The tempera- turewassampledrandomlyfromauniformdistributionwithmin- imumandmaximumvaluesof20^◦Cand60^◦Crespectively.Theun- certainty bands, 10th percentile, mean and 90th percentile were

used to quantify the degree of uncertainty for each of the pre- dictedmodelstatesresultingfromtemperatureuncertainty.These importantlevelswithintheuncertaintyrangewillbeusedtoprop- agatethestatesuncertaintyontotheboundaryoftheattainablere- gionsinordertodefinetheself-optimizingattainableregions(see Section4.3).

3.2.3. Attainableregionanalysis

After defining the process model and calculating the output (model states) uncertainty bands resulting from input (tempera- ture) uncertainty, attainable region analysis can now be used to definetheperformancetargetsoftheanaerobictreatmentprocess.

The following section outlines the methodologicalframework for ARconstruction,itsapplicationforprocesssynthesisaswellasfor definingtheself-optimizingoperatinglimitsundernon-isothermal conditions.Theframeworkinvolvesfivemainsteps:

Step1:Preparation

Thisinvolvesdefinitionofthereaction kinetics,AR dimension, statevariables(thoseusedtorepresenttheAR)aswellasthefeed point andtemperature valuesthat correspond to the mean,10th and90thpercentileofthestateuncertaintybands.Astoichiomet- ricschemeofthebioreactionoccurringinthemethanebioreactor consist oftwo main reactions catalyzedby acid-formingbacteria, Eq.(16)andmethane-formingbacteriaEq.(17)

k1SBVS rXac

→ Xac+k2SVFA (16) k3SVFA

rXme

→ Xme+k4CH4 (17) Letting rows 1–5 correspond to S_BVS, X_ac,S_VFA,X_meandCH₄ re- spectively, the stoichiometric coefficient matrix A is therefore a 5×2matrix,givenbyEq.(18)

A=

⎡

⎢ ⎢

⎣

−k1 0

1 0

k₂ −k3

0 1

0 k4

⎤

⎥ ⎥

⎦

⁽¹⁸⁾

Sincethere aretwo independentreactions participatinginthe system(Rank(^A)=2),weexpectthesetofpointsgeneratedbythe anaerobictreatment process to reside ina two-dimensional sub- space in R⁵. As all model outputs are functions ofvolatile fatty acidsandconcentrationofmethanogenicbacteria,itissensibleto generatetheARin(S_VFA−Xme)space,whichprovidesinformation requiredtomaximizegasproductionandvolatilesolidsreduction.

The number of dimensions in which the AR must be con- structedwasreducedusingtheconceptofyieldcoefficients,which hasbeenusedpreviouslytoreducethenumberofdimensionsdur- ingARanalysis(Scottetal.,2013).

ThisimpliesthattheconcentrationsofBVSandacidogenicbac- teriacanbeexpressedasafunctionofVFAandmethanogenicbac- teriaconcentrationsasinEqs.(19)and(201).

Xac=Xac0+ 1 k2

SVFA−SVFA0+k3

(

^Xme−Xme0

)

(19)

SBVS=SBVS0−k1

k2

SVFA−SVFA0+k3

(

^Xme−Xme0

)

(20) The ability to calculateXac andS_BVS asa function of Xme and S_VFAallowustoalsoexpresstherateandconcentrationvectorsof

XmeandS_VFA exclusively.Inother words,foreach Xme andS_VFAin the(S_VFA−Xme)spacewecancalculatearatevectorthatuniquely determinestheCSTR locusandPFRtrajectoryfromaspecifiedor- ganicload.

Step2:ARconstruction

Foreachofthethreetemperaturevalues(thetemperatureval- uesthat correspondtothemean,10thand90thpercentileofthe state uncertainty bands), the AR is generated using a combina- tionofPFR,CSTR andmixing. Thisisthemostdifficultandtime- consuming step butalso provides the mostvaluable information abouttheoperatinglimitsofthesystem.ARconstructiontypically beginsbydeterminingthePFRtrajectoryandCSTRlocusfromthe feed.The PFR trajectory isthe setof pointsgeneratedby solving the steadystate model ofa PFR reactor(aset ofordinary differ- ential equations) whiletheCSTR locusisthe setofpointsgener- atedbysolvingtheCSTRmodel(asetofnonlinearequations).The convexhullforthesetofallpossiblepointsgeneratedbyallpos- siblecombinationsofPFR,CSTRandmixingdefinestheattainable region.The CSTRequationsare solvedusingNewton method,im- plemented bythe Matlabroutine ‘fsolve’whilethePFRequations are solved usingthe Runge-Kutta4th to5thorder algorithm im- plemented by the Matlabode45 routine for solving non-stiff dif- ferential equations.The convexhalloftheentiresetofgeometric pointsisobtainedbyusingtheMatlab‘convhull’routine,whichim- plementstheQhullalgorithm(Mathworks,NatickNA).

Step3:BoundaryInterpretation

ThisstepinvolvesinterpretationoftheAR boundaryintoreac- torstructures, basedon thefundamentalcharacteristicsoftheAR boundary. The boundary of the AR is composed of reaction and mixingsurfacesonly.Reaction surfacesarealwaysconvexandthe pointsthatformconvexsectionsoftheARboundaryarisefromef- fluentconcentrationsspecificallyfromPFRtrajectories.Foratwo- dimensionalsystem,pointsontheARboundarythatinitiatethese convex PFR trajectories arise from specialized CSTRs. This infor- mation is used to determine digester configurations required to achieve thepointslocatedwithinandontheboundary oftheat- tainableregion.

Step4:Defineself-optimizingoperatinglimits

Afterconstructingtheattainableregionsthatcorrespondtothe mean,10thand90thpercentileofthestateuncertaintybands,the self-optimizingoperatinglimitsofthesystemcanbeobtained.This isdonebyoverlayingeach ofthethreeARsontooneanotherand determiningtheregionofintersectionbetweenalltheregions.The intersectionregion, thoughusually smallerthan eachofthe indi- vidual ARs, willalways be attainable inspite ofthevariationsin temperaturewithinthepredefineduncertaintyrange.Sincetheen- tireboundaryoftheindividualARshavealreadybeeninterpreted in terms of reactor structures (step 3), the particular reactor re- quiredtoachievepointsontheself-optimizingARisknown.

4. Resultsanddiscussion 4.1. Temperaturesensitivityanalysis

As earlier pointed out in the introduction to this article, the goal hasbeen to define theperformance targetsofbiogas plants operated undernon-isothermal conditions(plants where temper- ature is not controlled). This showsa need to be explicit about exactly what ismeant by theword performance targetsfornon- isothermal biogas plants,which implieswhat isachievable albeit all possiblevariationsintemperature.Figs. 3and4comparesthe resultsobtainedfromanalysingthesensitivitiesofthefourprocess statestotheoperatingtemperatureofthebiogasdigester.

From the Figures, it can be seen that within the 28 days of anaerobic digestion, thegreatest influence oftemperatureon the processstatesisduringstart-up(firstfivedaysofoperation).How-

ever, at this stage, a conclusion cannot be made with certainty, aswe don’t knowif thiseffectcould be dueto definedvalue of temeprature(35^◦C)usedtoobtainthesensitivityfunctions.Recall fromSection3.1.2thatthetemperaturedependenceoftheanaero- bictreatmentprocesshasbeenmodelledusingtheRatkowskyex- pandedsquarerootmodel,Eq.(11).Becausethereisnoexplicitso- lutionforthekineticmodel(systemofordinarydifferentialequa- tions) of the anaerobic treatment process, the sensitivity equa- tions are coupled with the original model differential equations andsolved to obtain the temperaturesensitivities for the neces- sarytimepoints.Thisentaildefiningavaluefortemperatureinthe Ratkowskyexpandedsquarerootmodelbeforethecoupledsystem ofequations canbe solved. In orderto ensure the obtainedsen- sitivityeffects are not depended upon the definedvalue oftem- perature,theMonte Carlo procedurewasapplied tovisualize the effect of inputs (different defined temperatures between 20 and 60^◦C) on the outputs (behaviorof the sensitivity functions). Fig- urefour presentstheoverall behaviorsofthesensitivityfunction for1000MonteCarlosimulations.Noticethatdespitevariationsin thedefinedvalueofthetemperature,thesensitivityoftheprocess statestotemperatureisstillverysignificantonlyatthestartupof theanaerobictreatmentprocess.

Takentogether,theobservationsfromFigs.3to5suggestthat temperaturehasa significant effectduring start-up ofthe anaer- obic treatment process. Interestingly, this correlation is because duringthe early stagesof anaerobic digestion, the anaerobicmi- croorganismsarestilltryingtoacclimatizetotheconditionsofthe waste. The acclimatization of anaerobic microorganisms is influ- encedbytemperatureandotherfactorssuchaswastecharacteris- tics,inoculumactivity,pH,loadingrate,retentiontimeandreactor configuration(WeilandandRozzi, 1991; Ghangrekar etal., 1996). This studyfocuses on the interactions betweentemperature and reactorconfigurationas thisaccordwith observationsfromother studies, investigatingthe effect of modifying the digester config- urationsto overcome the influenceof temperature (Chong etal., 2012;Gabyetal.,2017).

4.2.Uncertaintyquantification

Havingdiscussedhow sensitivethe process statesare totem- perature, we now present the actual quantification of the over- all effect of temperature on the process states. Fig. 6 presents thepropagationoftemperatureuncertaintyontothemodelstates, whichhasbeenquantifiedusingthe10thpercentile,meanaswell as the 90th percentile. Observe that the uncertainty band (10th and90thpercentile) is widestwithin thefirst few days ofoper- atingtheanaerobictreatmentprocess,whichfurthersupportsthe high temperature influence during process start-up. It is impor- tant for readersto note that the uncertainty bands presented in Fig.6 are dependenton the temperatureuncertaintyrange (20 - 60^◦C)usedfortheMonteCarlosimulation.

The results are interpreted as follows: the larger the uncer- tainty band (difference between the 10th and 90th percentiles), thelower themodel prediction quality.Accordingly, dueto tem- perature uncertainty, the model prediction quality is low for all thefourstateswithin thefirst fewdays ofoperation.The attain- able regions, which defines the limits of achievability of a pro- cessisbasedon theattainabilityofstates anduncertainty inthe stateswillthereforeresultsinuncertaintyindefining theoperat- inglimitsofthesystem. Hence modeluncertainty,resultingfrom temperatureuncertaintyreducesthereliabilityofusingattainable regions for defining process performance targets. This is where the strength of this study comes into play where we quantify themodel predictionuncertaintyandincorporateits effect when definingtheoperatinglimitsofthesystem(objectofthenextsec- tion).

Fig. 3. Sensitivity functions of model states to temperature.

Fig. 4. Superposing sensitivity functions of model states to temperature.

4.3.Operatinglimitsandoptimalplantconfigurations

So far, this paper has illustrated the sensitivity of the model statestotemperatureandhowuncertaintyintemperaturereduces

thequalityofmodelprediction.The followingsection willdiscuss how to correct formodel prediction uncertaintywhen using the model to define the operating limitsof the anaerobic treatment process.Theattainableregionisdefinedastheconvexhullforthe set ofstatesachievable fromdifferentcombinations ofPFR, CSTR andmixing. Hence inorder to factorin temperatureuncertainty, wemakeuseofthestate’suncertaintybands(mean,10thand90th percentile)thathavebeencomputedinSection4.2.Thisisdoneby constructing the attainable regions using the temperature values that correspondto uncertaintybandsanddeterminingthe region ofintersectionbetweenthethree individual regions.The temper- ature valuesthe gave the mean,10th and90thpercentile of the statepredictionwasrespectively31.60,20.00and52.40^◦C.

Figs. 7 to 9 present the attainable regions constructed using thetemperaturesthatcorrespondtothe10thpercentile,meanand 90th percentile state uncertainty bands respectively. The bound- ariesoftheattainableregions havebeeninterpretedintodigester configurations,whichcanbeusedtoattainthedifferentoperating limitsdefinedby the regions. The boundaryof theattainable re- gionrepresentsthesmallestsubsetofpointsthatcangenerateall otherpointsachievablebythesystemusingpossiblecombinations offundamentalreactortypesandmixing(Ming etal.,2016).This smallestsubset ofpointsis referredto astheconvexhull, which fora two-dimensional AR,is interpreted as thesmallest polygon enclosed by planar facetswhereby all elementslie onor within thepolygon(Asieduetal.,2015).Theauthorsnowillustratethein-

Fig. 5. Influence of temperature on process states for 10 0 0 Monte Carlo Simulations.

terpretationoftheARboundaryintoreactorstructuresusingFig.7. OnFig.7,thepointAisthefeed,whilethelineACandregionde- finedbypointsABCistheAR.TheconvexsegmentABisknownas thePFRtrajectory whilesegment AtoDiscalledCSTR locus.The curvesrepresentedby E(movingfromCSTR locusto pointC)are PFR trajectories obtainedusing concentrations onthe CSTR locus asfeedconditions.ThepointCisthereforeobtainedbyrunninga CSTR fromthefeed(pointA) followedby aPFR fromCSTR efflu- ent. As mentioned in Section 3.2.3 (step 3), the boundary ofthe AR iscomposed ofreactionandmixingsurfacesonly.Mixingsur- faces are always straight lines while reactionsurfaces are always convexandthepointsthatformconvexsectionsoftheARbound- aryarisefromeffluentconcentrationsspecificallyfromPFRtrajec- tories. The lineACis thereforea mixing linewhilethe curve AB isareactionsurface.Concentrationsthatliealongthemixingline AC(C_AC) canbeobtainedbymixingpointsAandC,andgenerally followsthelever-armrule,Eq.(21).Thereactorstructurerequired to achieve points on the mixingline AC is therefore given by a CSTR+PFR(pointC)withabypassfrompointA.

CAC=

α

^CA+

(

¹⁻

α )

^{CC, 0≤}

α

^{≤1 (21)}

Where

α

^{isknownasthemixingratio.Similarreactorinterpre-}

tationsweremadefortheother substratesaspresentedinFigs.8 and9.Itisnottheintentionofthisarticletogointodetailedex- planation ofthegeometryinvolvedininterpreting theARbound- aryintothedigesterconfigurations,interestedreaderscanconsult thepreviouslycitedpapers.

Three interesting remarks can be madefrom theFigs. 7 to 9. (1)Anincreaseintemperatureofthedigester increasestheoper- ating limitsof the system. The operating limit is definedby the area of the convex hull, which is computed as of 1.53×10⁻⁴, 4.95×10⁻⁴ and6.32×10⁻⁴ (^g/L)^{2 forthe 10th}percentile, mean and90thpercentilerespectively.Thisbecauseanincreaseintem- perature(withinacertainrange)generallyincreasestherateofthe anerobic digestion and hence more states will be output by the systemsoperating athigher temperatures.(2) Usinga multistage digesterconfigurationasopposedtoasingledigesterincreasesthe operatinglimitsoftheanerobictreatmentprocess.Pointsthatcan be achievedby the system are onlythose located within the at- tainable region.Forall thetemperatureshigher concentrationsof methanogenic archaea can be obtainedfor higher concentrations ofvolatile fattyacids only using a multistage digester configura- tion consisting of a CSTR followed by a PFR and by pass from feed.This is because multistage digester configurations are opti- mizedeachstep(acidformationandmethaneproductionsteps)of theanaerobictreatment process(EPA,2006). (3)Achange inthe operatingtemperaturedoesnotaffectthegeometryoftheattain- ableregion boundary butonly thelimits ofachievability, area of the convexhull (as explained in the first remark above). This is quiteinteresting becausethegeometryoftheattainableregion is uniqueforagivenkineticsandfeedpointandnotfortemperature (HildebrandtandGlasser,1990;Hildebrandtetal.,1990).

Theresultsallputtogetherimplythattemperatureanddigester configuration have a significant effect on the operating limits of

Fig. 6. State uncertainty bands obtained from 10 0 0 Monte-Carlo Simulations.

Fig. 7. Attainable region (operating limits) for 10th percentile state prediction.

the anaerobic treatment process and hence modifying the tem- peratureorreactorconfigurationwouldaffecttheperformance of theanaerobictreatmentprocess. Thisaccordsthefindings ofpre- vious studies reporting that modifying the reactor configuration andor temperature influences the performance of the anaerobic treatmentprocess. Apractical exampleis theuseof amultistage dieter consisting of a UASB (a kind ofplug flow digester) anda

Fig. 8. Attainable region (operating limits) for mean state prediction.

CSTRtoimproveperformancetargetsofanaerobicdigestionunder variable temperatureconditions(Lettinga andHulshoff Pol,1991; Mahmoudetal.,2004).

Fig. 10present the intersectionof the threeregions to define the self-optimizing attainable region of the anaerobic treatment process. Notice that the intersection region (region bounded by

Fig. 9. Attainable region (operating limits) for 90th percentile state prediction.

Fig. 10. Superposition of attainable regions showing the self-optimizing attainable regions.

ABC)corresponds to the10thpercentile attainable region.Thisis an interesting featurethat theauthors noticewhendefining self- optimizing attainable regions using temperature uncertainty (as opposed to kinetic uncertainty). An explanation forthis observa- tion is that the statesthat are achievable at lower temperatures willalsobeachievedathighertemperaturesbutnotallstatesthat are achievable at higher temperatures can be achieved at lower temperatures. Put it in another way, since the rate of digestion increases with temperature (of course within a given range), all other things beingequal, a givenquantity ofbiogas produced by a plant operating at a lower temperature can also be produced by a plantoperating ahigher temperature. Thedifference is that theplantoperatingathighertemperaturewillproducemorebio- gasinadditiontotheonethatisproducedbytheplantoperating at lower temperature. However, generally, the self-optimizing at- tainableregionisgenerallysmallerthanthetrueARbecauseself- optimalityimpliesan acceptablelossinoperatinglimits. Itisnec- essaryheretore-define exactlywhat ismeant byself-optimizing attainable regions. Unlikethe attainable region, whichrepresents the setof all possiblestates that isachievable by thesystem for a definedtemperature,kineticsandorganicload(feedpoint), the self-optimizingattainable regionrepresentstheset ofallpossible states attainableby the systemevenin casesof temperatureun- certainty.The sizeoftheself-optimizing attainableregionis con- ditionaltotheuncertaintydomain(20– 60^◦C)definedfortemper-

atureduringthepropagationoftemperatureuncertaintyusingthe MonteCarlosimulationprocedure.Thelargerthedomainofuncer- tainty,thesmallertheareaoftheself-optimizingARandviceversa andhence incorporatingtemperatureuncertaintyreduces theop- eratinglimitsofsystem.However,otherthandefiningaverylarge operating limit, which will be achieved some of the times (for fixedtemperaturevalues),wedefineasmallerlimitwhichwillbe achievedallofthetimes(despitevariationsintemperature).

Thefindings fromthisstudyare thereforehighlyimportantin making economic feasibility decisions about the performance of biogasplantsespeciallyincaseswhereaccuracyisverynecessary.

Thisisakeymotivationforrenewableenergyinvestorsasitisbet- ter tomaking investmentdecisionsforbiogas plantsusingworse caseperformancetarget.Thisisbecauseifevaluationsindicatethat theinvestmentis profitable,then evenbetter profitability indices will be obtainedduring real-timeoperation if the environmental conditions(mainlytemperature)favorthekineticsoftheanaerobic treatmentprocess. Evenin caseswheretheenvironmentalcondi- tions are unfavorable, we are sure to still be profitable since in- vestment decisionshave beenmadeusing aworse-case scenario.

Theresultsofthisstudyarethereforehighapplicabletocountries wherebiogas plantsareoperatedeitherunderisothermalornon- isothermalconditions.

Summarily,the results, which show that the operating limits, definedbytheARdiffersforeachoperatingtemperature,provides thefollowingnoteworthycontributions ofpracticalrelevance:For plants operated under isothermal conditions, the study provides aframeworkforobtainingunique optimaldigester configurations, which are specific to the digestion temperature under consider- ation (psychrophilic, mesophilic or thermophilic ranges). For the caseofnon-isothermalplants,thestudyproposesaself-optimizing approachtodefine performancetargets, whichare optimalforall temperatures within a defined boundary. This studythough pre- liminarypresentsabreakthroughinextendingtheuseofdigester networkstosolvemoreoperationalchallengesassuchsystemscan optimizeeverystepintheanaerobictreatmentprocess.Foranal- readyexistinganaerobicdigestionplant,theARconceptshowsthe proximity of the existing systemin relation to the absolute best performance,whichisimportantindecidingwhetherornottoin- vestresourcestooptimizetheplant.

It is interesting to compare the approach to digester network synthesispresented inthisstudywiththat presentedinprevious studies. The conventional superstructure optimization technique fordigesternetworksynthesisPontesandPinto,2009)suffersfrom multiple solutions (or local optimum) because it involves a very largereactor superstructureandprovides nosystematic approach foranswering thefollowingthreequestionsabouttheoptimaldi- gesternetwork:((1)whatnumberofindividualdigestersshouldbe inanoptimalnetwork.(2)Shouldby-passorrecyclestreamsbein- cluded,ifyes(3)whereinthenetworkshouldtheybepositioned.

TheuniquenessandstrengthoftheARapproachisthatitprovides theabsolutebestperformance(totalityofachievablestates)forall possibledigesterconfigurations(eventhosethathavenotyetbeen devised).

As mentioned in the introduction, the other papers on the AR concept work for isothermal biogas plants. This study is highly novel by introducing a self-optimizing concept for non- isothermalconditionsthroughthepropagationoftemperatureun- certainty onto the AR boundary to define the self-optimizing at- tainable regions. The idea of self-optimizing attainable regions first introduced in a recent paper published by the authors (Abunde Neba et al., 2020) is that instead of having an optimal performancetarget,whichcanonlybeachievedsomeofthetimes (duetotemperaturefluctuations),itisbettertodefineanearopti- maloperatingtarget,whichisachievableallofthetimes.Thecur- rentpaperhasextendedtheconcept ofself-optimizingattainable