Biodigester rapid analysis and design system (B-RADeS): A candidate attainable region-based simulator for the synthesis of biogas reactor structures

Anaerobic digesters are seldom designed based on process kinetics, but rather on a combination of hydraulic and organic loading, which may limit operational performance. This study focuses on the incorporation of process kinetics in the design of anaerobic digesters, within the attainable region conceptual framework. Candidate attainable regions for anaerobic digesters are identified using the software environment Biodigester Rapid Analysis and Design System (B-RADeS), which couples, biodegradation kinetics as well as economic parameters for the synthesis of biodigester structures. By considering swine, palm oil and pharmaceutical wastewaters, payback periods of 0.5, 1 & 2 years, and substrate, kinetic model and/or economic parameters, a promising digester structure (and associated hydraulic retention times) is synthesized, consisting of a CSTR followed by PFR (15 days), CSTR (4.8 hours) and a PFR with bypass of feed (3 days). The framework offers great promise for widespread practical application. © 2019 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license. ( http://creativecommons.org/licenses/by/4.0/ )


Introduction
Anaerobic digestion of solid waste and/or wastewater sludges has long been used for stabilization of these wastes prior to disposal. Among the benefits involved in anaerobic waste treatment compared to the aerobic counterpart are: improved dewaterability of the treated waste and generation of renewable bioenergy ( Mes et al., 2003 ). The construction and operation of anaerobic treatment plants requires optimizing the techno-economic feasibility by defining optimal process configuration of anaerobic digesters. There exist several types of anaerobic digesters each of which have specific characteristics making them more adequate to treat specific types of organic wastes ( Mao et al., 2015 ). For treatment of solid waste and sludges, low-rate anaerobic systems are more appropriate due to their use of long but coupled hydraulic and sludge retention times to ensure a stable operation of the pro-

B t
Annual savings from electricity consumption ($) C Gen Cost of biogas electricity generator ($) C Inv Investment cost ($) C con Cost of digester construction ($) C m Annual cost of digester maintenance ($) C misc Cost associated to miscellaneous activities ($) C pf Annual cost of biogas purification ($) C t Annual operating cost ($) K * Kinetic constant of Chen and Hashimoto based on specific methane production rate (−) K ha Maximum rate of substrate degradation by the acidogenic bacteria ( da y −1 ) K Lr Monod-like half saturation constant for continuous mode operation ( gVS / L / day ) K am Maximum rate of substrate utilization by the acetogenic/methanogenic microorganisms K i Inhibition constant for biogas yield ( g / L ) K s Monod-like half saturation constant for continuous mode operation ( gVS / L ) L r Organic loading rate ( gVS / L / day ) P el Percentage of methane utilized for electricity generation Pr pf Price for biogas purification per unit volume ($/ m 3 ) P t Net annual benefit ($) Recalcitrant fraction of initial volatile substrate that is non-biodegradable.

R p
Specific rate of biogas production by acetoclastic methanogens ( mL biogas / gVS / day ) R pm Maximum specific rate of biogas production by acetoclastic methanogens ( mL b / gVS / day ) S i Initial concentration of substrate taken up by acetogenic/methnogenic micororgansims ( g / L ) S o Initial substrate concentration ( gVS / L ) S u Concentration of acidified substrate produced by acidogenic bacteria ( g / L ) T el Annual time period for use of electricity ( d ) T pr Annual time period for biogas purification ( d ) V D Volume of digester (mL) X T X Characteristic matrix X am concentration of acetogenic/methanogenic microorganisms Y PS Biogas yield coefficient ( Rate of change in biogas yield −d S/d t Rate of decrease in the concentration of hydrolyzed substrate ( gVS / L / day ) k i Inhibition constant for cell growth( g / L ) k n ( s ) Rate of substrate degradation by acidogenic bacteria ( g / g / day ) k s Monods' half saturaion contant for substrate uptake m VS Mass of volatile solids added into the digester (gVS) n t Project lifespan (years) r p Modified rate of biogas production by acetogenic/ methanogenic microorganisms ( mL b / gVS / day )  sludges and solid wastes. Using empirical methods to optimize design of anaerobic digesters often requires construction of expensive prototype systems and time consuming studies, which has been a key motivation for reliance on model-based techniques ( Yu et al., 2013 ). Use of the models is again highly dependent on the availability of kinetic coefficients and hence modelling requirements for design of anaerobic digesters are often simplified to a minimal number of inputs and experimental states (most commonly biogas yield) required for model identification ( Batstone, 2006 ). The simplified models employed for digester design can be single stage, based on the rate limiting step approach, or two-stage, based on lumping the process in to acid-forming and methaneproducing microorganism. The single stage models previously employed in digester design include first order models ( Linke, 2006 ;Momoh and Nwaogazie, 2011 ) and the models based on maximum bacteria growth rate ( Fdez.-Güelfo et al., 2011 ;Fdez-Güelfo et al., 2012 ) while the two stage models include the biogas yield models presented by Momoh et al. (2013) . Although simplified models have been used extensively in digester design, published articles are limited to mainly to determination of digester capacity based on parameters such a VS loading, temperature, etc. with first order models being mostly used ( Momoh et al., 2013 ;Wang et al., 2007 ). For design, critical issues are hydrodynamics, as well as the behaviour of solids, which requires at least two-stage, with hydrolysis and biological steps ( Batstone, 2006 ). This study focuses on hydrodynamics and an approach to optimize hydrodynamic configuration of anaerobic digesters based on simplified models will be a highly practicable as it will require less experimental measurements to estimate kinetic constants. Although the study develops two-stage models which can be identified based using only biogas yield measurements, the emphasis of the paper is not necessarily on the models but on how the authors use the models to develop new hydrodynamic configurations for operating low rate anaerobic digesters. The design objective is to minimize the process time as well as payback period by considering biodegradation and mixing as the only permitted fundamental processes occurring in the digester. The approach is based on the concept of attainable regions (AR), which is a technique for process synthesis and optimization that incorporates elements of geometry to understand how networks of chemical reactors can be designed and improved. "Following this initial work, many other researchers advanced AR research. Glasser et al. (1987) proposed a geometric approach that identified candidate AR's satisfying a number of necessary conditions that the AR must possess. Burri et al. (2002) demonstrated that, within the Infinite DimEnsionAl State-space (IDEAS) conceptual framework, construction of the true AR, and increasingly accurate AR approximants, can be carried out through Infinite Linear Programming (ILP), and a sequence of approximating finite Linear Programs (LP) respectively. Subsequently, Manousiouthakis et al. (2004) developed, within IDEAS, necessary and sufficient conditions that the true AR must satisfy, proposed the Shrink-Wrap algorithm for AR construction, and established, this algorithm's equivalence to the aforementioned LP based AR construction methods. They also demonstrated that the true AR can be potentially larger than the candidate AR's identified by geometric methods. Zhou and Manousiouthakis (2006) demonstrated that the true AR for reactor networks involving only reaction and mixing may be smaller than the true AR for reactor networks also incorporating diffusion effects (e.g. by considering non-ideal dispersion reactor models). Zhou and Manousiouthakis carried out pollution prevention studies using the AR approach ( Zhou and Manousiouthakis, 2007a ), extended the AR approach to reactor networks involving variable density fluids ( Zhou and Manousiouthakis, 2007b ), discussed the dimensionality of the space in which AR construction can be pursued ( Zhou and Manousiouthakis, 2008 ), and extended the AR approach to nonisothermal reactor networks ( Zhou and Manousiouthakis, 2009 ). Around the same time, Posada and Manousiouthakis (2008) , proposed AR construction methods for reactor networks with multiple feeds, while Davis et al. extended the AR approach to batch reactor networks ( Davis et al., 2008 ). More recently, Conner and Manousiouthakis extended the AR approach to general process networks ( Conner and Manousiouthakis, 2014 ), while Ming et al. (2016) summarized many of the AR literature results." Geometrically, the attainable region represents the region bounded by the convex hull for the set of points achievable by the fundamental processes occurring in the system ( Asiedu et al., 2015 ;. Once the AR has been determined, the limits of achievability by the system for the given kinetics and feed point is known and the boundary of the AR can then be used to answer different design and/or optimization questions related to the system. Our recent publication, Abunde Neba et al. (2019) has been first of its kind laying down theoretical framework for use of attainable regions to model optimal configurations of multistage anaerobic digesters. The study employed a four-state dynamic model of anaerobic treatment process and the attainable region analysis has been based on concentration (state) space but not residence time. The lack of residence time makes it impossible to size the digester structure or perform economic feasibility studies on the optimal digester structure. In addition, the four-state model (compared to the simplified model in this study) poses requirement for more experimental measurements hence limiting its application to situations where process measurements are limited. The current study is designed to illustrate how simplified models (requiring only biogas yield measurements) of the anaerobic treatment process can be used for attainable region analysis involving residence time space. "AlHusseini and Manousiouthakis (2013) were the first to incorporate residence-time considerations in the AR conceptual framework, by first introducing a production normalized, capital cost measure for a reactor network, that they termed "Network Residence Time" (NRT), and defined as "the ratio of the sum of the volumes of all reactors participating in the reactor network over the total volumetric flowrate entering the network." They subsequently introduced the Network Residence Time Constrained Attainable Region (NRT-C-AR), which they then proceeded to quantify using a Linear Programming formulation within the IDEAS conceptual framework. The advantage of constructing an AR of the anaerobic treatment process that incorporates residence time considerations, over the AR presented in our previous study, is that it enables the coupling of biodegradation kinetics, economic feasibility objectives and country specific macroeconomic parameters for the synthesis of biogas digester structures. By its use of attainable regions, knowledge of all possible states, for all possible digester configurations can be obtained considering biodegradation and mixing as the only fundamental processes occurring in the digester. Unlike previous studies where economic analysis is performed to determine the feasibility parameters of a predefined digester configuration, this study rather determines the biodigester network configuration required to achieve a given economic objective based on the macroeconomic situation of a given country. Finally, the study seeks to deploy the theoretical framework into a software in order to save time and effort for designers who are planning and designing biogas plants for different process or economic scenarios.

Attainable region theory for process synthesis and optimization
The Attainable Region (AR) theory is a technique that incorporates elements of geometry and mathematical optimization, to design and improve operation of chemical reactors ( Ming et al., 2016 ). The power of the AR approach to process optimization is that the answer to all possible optimization problems, even the ones not considered are first determine, and then we look for ways of achieving that answer. In reactor operation knowledge of all possible reactor states for all possible reactor configurations, even those that have not yet been devised, is obtained. For a two-dimensional system, the convex hull for the set of all points achievable by all possible combinations of CSTR + PFR and mixing defines the attain-able region. For higher dimensional systems, the attainable region is the convex hull for the set of points generated by all possible combinations of CSTR, PFR, DSR and mixing lines. The convex hull is understood as the smallest subset of a set of points that can be used to generate all other points by reaction and mixing ( Ming et al., 2016 ). Geometrically, a convex hull is a finite convex polytope enclosed by a finite number of hyperplanes, which is interpreted in a two-dimensional space as the smallest polygon enclosed by planar facets such that all of the elements lie on or in the interior of the polygon ( Asiedu et al., 2015 ). Once the AR has been determined, the limits of achievability by the system for the given kinetics and feed point is known, which can then be used to answer different design or optimization questions related to the system.
Given a set of reactions and associated kinetics, the following five key steps needs to be performed in order to complete an attainable region analysis ( Ming et al., 2016 ): ➢ Define the reaction, dimension and feed set. ➢ Generate the AR using combinations of the fundamental processes. ➢ Interpret the AR boundary in terms of reactor equipment. ➢ Define the objective function and overlay this onto the AR to determine point of intersection with the AR boundary. ➢ Determine the specific reactor configuration required to achieve the intersection point.
Some necessary conditions for AR derived from the work of Glasser et al. (1987) can be summarized as follows: ➢ The AR includes all feeds to the system. ➢ The AR is convex. ➢ No process vector point out of the AR boundary. ➢ No rate vectors in the complement of the AR when extended backward intersects the AR.
The objective of this section is to analyze the aforementioned necessary requirements with respect to its application to the anaerobic treatment process. However, AR analysis requires that the process kinetics is known and we therefore begin by modeling the kinetics of the anaerobic treatment process.

Reaction kinetics of the anaerobic treatment process
In the present paper, the mathematical models describing the kinetics of substrate utilization and methane production in anaerobic treatment process are developed based on the approach presented by Momoh et al. (2013) . The approach assumes that the AD process takes place in three stages. (i) hydrolysis/acidogenesis of the organic substrates in wastewater by acidogenic bacteria to produce acidified substrate; (ii) uptake of acidified substrate by acetogenic/methanogenic bacteria and (iii) acidified substrate assimilation, growth and biogas production by the acetogenic/methanogenic bacteria. Fig. 2 presents the algorithm used to develop and validate the simplified two-stage based modes to predict biogas yield. The model development involves five main aspects, which include:

Enunciation of the process model
➢ Develop substrate degradation model. ➢ Formulate substrate uptake model. ➢ Choose microbial kinetic model. ➢ Derive models for substrate assimilation and biogas production. ➢ Identification of the developed model.
Considering these aspects led to a series of ordinary differential equations to predict biogas yield based on microbial growth kinetics Stage 1: Hydrolysis and acidogenesis Many constituents of organic wastes behave as complex substrates (polysaccharides, proteins, fats etc.). The Grau model presented in Eq. (1) , which has widely been used to model multiple substrate removal kinetics ( Kim et al., 2006 ;Liu, 2006 ) was therefore adopted for this study.
where −d S/d t represents the rate of decrease in concentration of substrate being hydrolyzed, S is the concentration of initial substrate left at every instant following onset of hydrolysis, S o is the concentration of initial substrate, k n ( s ) is the rate of substrate degradation by acidogenic bacteria, X is the concentration of acidogenic bacteria and n defines the degree of adaptation by acidogenic bacteria for substrate degradation.
The multicomponent substrate degradation model is based on the assumption that the different components are simultaneously removed and transported into the cells ( Grau et al., 1975 ). Assuming hydrolysis and acidogenesis are catalyzed by acidogenic bacteria, whose concentration is constant, then Eq. (1) can be re-written as Eq. (2) .
Where, K ha is the maximum rate of substrate degradation by acidogenic bacteria. Since anaerobic digestion is a biological process and the AR approach considers biodegradation and mixing as the only fundamental processes occurring in the digester, it becomes important to consider the non-biodegradable part of the substrate. The model is then modified as shown in Eq. ( Eq. (3) represents the kinetics of substrate degradation, where R f is the recalcitrant fraction of initial volatile substrate that is nonbiodegradable.
This model takes into consideration the acidified substrate produced after substrate degradation by acidogenic bacteria as well as the uptake of acidified substrate by acetogenic/methanogenic microorganism. S u represents the actual amount of the substrate that was acidified and utilized by the acetogenic/methanogenic bacteria while b is the fraction of initial volatile solids remaining in ef- represents the rate limiting coefficient for very slow (case of 0 <α< 1) or very fast (case of α = 1) metabolism of acidified substrate by the acetogenic/methanogenic bacteria ( Momoh et al., 2013 ). The constant K am is the maximum rate of substrate utilization by the acetogenic/methanogenic microorganisms.
Stage 3: Kinetics of bacteria growth and biogas production The attainable region is unique for a given kinetics and a change in organic substrate can cause a change the kinetic model used to describe the growth of microorganisms. Table 1 presents a list of microbial growth models considered to model substrate assimilation. The table has been assembled from Kythreotou et al. (2014) , who presented a comprehensive review of simple to scientific models for anaerobic digestion. As expected, the different models have different characteristics often making Integrates effect of microbial adoption to stationary processes by mutation Tessier model An exponential function used to describe cell growth processes For growth process affected by the allosteric effectors present in the acidified substrate Haldane model for product inhibition Moser, 1981und Bergter, 1983 Production inhibition model with effect of microbial adoption them more adequate to describe microbial growth of specific effluent and/or digester conditions rather than others. k s is the Monods' half saturaion contant for substrate uptake, μ max is the maximum specific growth rate for methatnogenic archae, m is the coefficient of acetogenic/methanogenic microbial adaptation for cooperativity, S i is the initial concentration of substrate taken up by acetogenic/methnogenic micororgansims, k is the kinetic constant of Chen and Hasshimoto, and k i is the substrate concentration where bacteria growth is reduced to 50% of the maximum specific growth rate due to substrate inhibition Taking the case of the Monod model for growth of acetogenic/methanogenic microorganisms, and using product and cell yield coefficients, the rate of biogas production can be expressed by Eq. (15) Where dy t / dt is the rate of change in biogas yield, y t is the biogas yield, Y PS is the biogas yield coefficient, Y XS is the cell yield coefficient and X am is the concentration of acetogenic/methanogenic microorganisms. If we consider the growth rate of the acetogenic/methanogenic bacteria is very slow or relatively constant while dy t / dt can be described as the specific biogas yield rate ( R p ) at the end of biogas production ( Momoh et al., 2013 ), then Eq. (15) can re-written as Eq. (16) . The parameter R pm = ( X am Y PS / Y XS ) is the maximum specific rate of biogas production.
Eq. (16) describes the biogas yield in anaerobic digester considering a batch operation mode. In cases where the system in operated in continuous mode, the initial substrate concentration ( S o ) is converted to loading rate by multiplying the factor ( Q / V ) as shown Table 2 Two-stage based models to predict biogas yield rate.
The resulting continuous mode counterpart of the biogas yield model considering Monod kinetics is shown by Eq. (18) .
Where, the variable L r is the organic loading rate into the biodigester and K Lr is the Monod's half saturation constant defined in terms of the organic loading rate. Similar process was applied to develop the other biogas yield rate models by assuming that the growth process of the acetogenic/methanogenic microorganism can be described using the other growth models presented in Table 1 . However, a parameter of k p was introduced to the Tessier based model as a coefficient to the exponential term, which serves as an index of the processing speed of R p as it approaches R pm due to the change in S o or L r (depending on the mode of operation).The derived biogas yield models considering a two-stage biodegradation kinetics is presented in Table 2 .
K is the kinetic constant of Chen and Hasshimoto defined in terms of specific biogas yield, K i is the substrate concentration where specific biogas yield rate is reduced to 50% of the maximum specific biogas yield rate due to substrate inhibition and n provides a useful measure of microbial cooperativity to biogas production.

Parameter estimation and statistical methods
The kinetic constants of the different models were estimated using the Matlab nonlinear regression solver ' nlinfit' (Mathworks Natick NA). In assessing the variability of the model identification process, we used the kernel density estimates and the parameter confidence regions. It is also interesting to note that marginal confidence intervals often used by several researchers do not account for correlations between the parameter estimates. Therefore, their use in parameter estimation can sometimes be misleading if there is strong correlation between several parameter estimates. In this study, we rather illustrate the use of joint confidence regions in assessing reliability of parameter estimates in least square regression.
The 100( 1 − α)% joint confidence region and the marginal confidence intervals of the parameter estimates is computed using Eqs. (29) and (30) respectively is the approximate standard error of the parameter estimates given by Eq. (31) is the deviation between the real ( β) and the estimated model paramters ( ˆ β ) , t υ,α/ 2 is the student's t-distribution parameter, ν is the number of degrees of freedom ( n − 1 ), where n is the number of data points used to compute the variance and average, F ( 1 −α) ,p, ( n −p ) is the F-distribution varince comparism parameter, p is the number of model parameters being estimated, ( n − p ) is the model degrees of freedom, and σ 2 is the true variance, and α = 0 . 05 is the level of significance.

AR analysis of the anaerobic treatment process 2.2.1. Reaction dimension and process vectors
Before it is possible to construct the AR, the engineer must first determine the space wherein the AR must reside (by choosing unique species components in the system that will represent the AR). These are variables required to characterize the state of the system, in this case, an anaerobic treatment process and must be sufficient to describe the dynamics of the fundamental processes chosen to describe the system. These variables would include biogas yield ( y t ) and retention time ( τ ). A key criteria for selecting variables in AR is that they must obey the linear mixing law. The concept of "Network Residence Time" (NRT), as introduced by Al-Husseini and Manousiouthakis (2013) defines the residence time of a reactor network as the ratio of the sum of the volumes of all reactor units constituting the network to the total volumetric flowrate entering the network. Using this definition, it can be shown that the residence time of a network comprising two reactor units obey the linear mixing law, Eq. (32) . This implies the overall residence time of the system must lie in a straight line between the residence times of the individual reactors, τ 1 and τ 2 comprising the system.
Where τ mix is the overall residence time of the system comprising two individual reactors. The developed simplified kinetic models predict biogas yield ( y t ), which is given in terms of volume of biogas produced (mL) per gram of volatile solids added to the digester (gVS). y t = V g / m V S . Assume we have two digesters of known biogas yield, we can obtain the actual volume of biogas produced for digesters 1 and 2 by V g1 = y t1 m V S1 and V g2 = y t2 m V S2 respectively. Conservation of mass may be used to calculate the total biogas yield for both digesters. Conservation of mass ensures that the total mass of volatile solids in the mixture is equal to the sum of the individual masses contained in beakers 1 and 2, which is given by m V ST = m V S1 + m V S2 . Computing the biogas yield of the entire system is equivalent to determining the biogas yield for a mixture of digesters 1 and 2 since the density of the liquid phase of the digester can be assumed constant. The biogas yield of the mixture ( y tM ) is given by the ratio of the total volume of biogas produced to the total mass of volatile acids added as shown by Eq. (33) .
If we set α = m V S1 / m V ST then Eq. (33) can be written as Eq. (34) , which is similar to the linear mixing law. What this means practically is that if we mix the contents of the liquid phase of two digesters, each of which contains a given quantity of volatile solids added, then the total biogas yield of the mixture will lie in a straight line joining that of both digesters.
The process of combining the contents of two parallel digesters (or digester networks) of different volatile solids contents results in a linear mixing law measured in term of biogas yield. This implies biogas yield may be used in the construction of candidate ARs in a similar manner to that for concentration.
The biogas yield and the retention time grouped together form a vector called the characteristic vector; C = [ y t τ ] , whose dimension determines the dimension of the optimization problem. We therefore have a 2-D optimization problem with the objective of minimizing [ τ ], time parameter.
If we assume that as substrate is consumed rate of change of biogas yield is directly correlated with the quantity of biogas to the biogas yield y t , such that the driving force for gas production is disappearing when the biogas yield gradually approaches its maximum ( y tm ) then for the mass of volatile solids added to the digester. This is modelled by Eq. (35) Where, r p is the modified rate of biogas production by acetogenic/methanogenic microorganisms. The reaction rate vector is therefore given by − − → r(C) = [ r p r τ ]

Generate the AR using combinations of the fundamental processes
The attainable region (AR) represents the set of all possible states that can be achieved by a combination two fundamental processes, biodegradation and mixing in the case of the anaerobic treatment process. In AR theory, mixing is performed by a continuous stirred tank reactor (CSTR) while biodegradation (reaction) is performed by the plug flow reactor (PFR), since the operation of both reactors respectively mimics the two fundamental processes. At steady state operation, the general mathematical model of a CSTR and PFR are given respectively by Eqs. (36) and (37) respectively.
C is the two-dimensional state vector made of biogas yield and the residence time, Eq. (38) while r ( C ) is the reaction rate vector, which can be illustrated to be given by Eq. (39) .
During construction of AR, the PFR trajectory is the set points generated by numerically solving the PFR equation while the CSTR locus is the set of points generated by solving the CSTR equation. The convex hull for the set of all points achievable by all possible combinations of CSTR + PFR defines the attainable region. The convex hull is understood as the smallest subset of a set of points that can be used to generate all other points by reaction and mixing ( Ming et al., 2016 ). Geometrically, a convex hull is a finite convex polytope enclosed by a finite number of hyperplanes, which is interpreted in a two-dimensional space as the smallest polygon enclosed by planar facets such that all of the elements lie on or in the interior of the polygon ( Asiedu et al., 2015 ) The candidate attainable region was constructed with Matlab using the following five-steps Step 1: Determine PFR trajectory from feed Step 2: Determine the CSTR locus from feed Step 3: Determine PFR trajectory from each CSTR point Step 4: Construct the convex hull of the set of achievable points Step 5: Verify the obtained AR against the necessary conditions of AR and if any condition is not met return extend the AR by running a PFR from the point of disagreement The PFR equations are solved using the Matlab ode45 routine for solving non-stiff differential equations while the system of nonlinear CSTR equations were solved using 'fsolve' routine The convex hall of the entire set of geometric points is obtained by using the Matlab ' convhull' routine, which implements the Quickhull algorithm (Mathworks, Natick NA).
It is important to mention that even though this construction approach has been applied in a couple of studies ( Ming et al., 2013 ;Ming et al., 2016 ), the NRT-C-AR obtained is a candidate and not the true NRT-C-AR. For a true AR, the Infinite DimEnsionAl Statespace (IDEAS) conceptual framework is applied to obtain a general linear programming formulation for the construction of the true NRT-C-AR, as shown in ( AlHusseini and Manousiouthakis, 2013 ). However, the interest of the study is not necessarily on the method used for AR construction, but on how the concept of attainable regions can be applied to optimized operation of the anaerobic treatment process. Also, even if just a candidate AR is obtained, it can still be used for process synthesis and optimization only that the totality of outputs is not obtained.

Interpret the AR boundary in terms of reactor equipment
The AR boundary is composed of two types of geometries: mixing lines referred to as lineations and manifolds of PFR trajectories referred to as protrusions. The role of PFRs on the AR boundary is to generate the outer extremities whereas CSTRs and DSRs (in the case of higher dimensional constructions) are used as connectors to these PFRs ( Ming et al., 2016 ). This implies the AR boundary is defined in terms of reactor structures, and for two-dimensional constructions, the boundary is composed of combinations of the two fundamental reactor types and mixing lines. The PFR and CSTR each exhibit unique geometric interpretation and hence determining the reactor configurations that form the AR relates to interpreting the surfaces that form the AR boundary using its geometric properties.

Define the objective function and the optimal reactor structure
The AR, which defines the limits of achievability by the system for the given kinetics and feed point can be used to answer one or more design or optimization questions related to the system. Twodimensional ARs involving residence time are particularly important for understanding the minimum reactor volume required for a given output. Since the construction and operation of anaerobic digesters generally requires capital investment, it would be interesting to use the AR concept in determining the profitability of the plant. However, we need to develop a suitable objective function that incorporates biogas yield and residence time (or digester volume).
The economic evaluation considers that biogas generated from the anaerobic digester is utilized for electricity generation. The total annual income, benefit ( B t ) from installing the biomethane plant is determined by Eq. (40) , which is benefit due to savings from electricity consumption.
Where P el is the percentage of methane utilized for electricity generation (which is 100% in the case of the current study), T el is the annual time period for use of electricity, b e is the unit conversion coefficient, Pr el is the feed-in tariff rate for biogas based electricity, gVS is the gram mass of volatile solids fed in the digester.
The total annual expenses or operating cost ( C t ) is computed by Eq. (41) . The operating costs are assumed to be a function of two factors: the repair and maintenance costs, Eq. (41a) , which is taken to by 1% of the cost of construction (0.01 C con ) and the cost of biogas upgrading, which is a function of the biogas volume, Eq. (41b) .
Where, C m is the annual cost of digester maintenance, C pf is the annual cost of biogas purification, C con is the cost of digester construction, T pr is the annual time period for biogas purification and Pr pf is the price for biogas purification per unit volume.
The cost of investment/construction is computed using the rates of a commercial biogas company in Ghana, stating the cost of digester construction to be $300 per cubic meter ( Mohammed et al., 2017 ). This includes administrative, transport costs, consultancy fees and other logistic aspects. The final expression of the total annual cost of digester operation is given by Eq. (42) .
Where V D is the volume of digester. The annual profit ( P t ), is defined as the difference between annual benefit, B t , due to savings from electricity consumptioin and the annual operating costs, C t . This is expressed by Eq. (43) .
Substituting the expressions for B t and C t into Eq. (43) the expression for the annual profit can be written as in Eq. (44) Since the AR is constructed in residence time space, it is necessary to express the volume of digester ( V D ) in terms of residence time τ and volumetric flow rate Q . The expression for P t as a function of residence time becomes; P t = gV S × y t 0 . 9 P el × T el × b e × P r el − T pr × P r p f − 3 τ Q The economic evaluation of the digester investment is based on the payback period (PBP) ( Gittinger, 1986 ) and the decision rule is that one generally accepts projects that require shorter number of years to recover the investment. The payback period is given by the annual profits, generated over n years, needed to recover the total C In v = C con + C Gen + C misc (46a) Where C Inv is the cost of investment, r is the discount rate, C Gen is the cost of biogas generator, C misc is the miscellaneous cost and n t is the project lifespan. Eq. (46) can be rearranged to express y t as a function of τ , Eq. (47) which may be plotted over the AR boundary as contours (for different specified values of n ) to determine the point(s) of intersection with the AR boundary. These intersection points represent the optimal operating point (which can interpreted into an optimal reactor structure) required to achieve a specified payback period. Table 3 presents of summary of the parameter sets that are used to perform the economic evaluation of designing a constructing a methane plant.

Development of computational model
The design of the graphical user interface (GUI) was done using the Matlab GUIDE (Graphical User Interface Development Environment). This is done using icon-based programing using several objects such as push buttons, static texts, edit texts, pop-up menus and axes handles. GUIDE generates a GUI and the m-file that contains the code to handle the initialization and launching of the GUI. After creation of the GUI, it was programmed by entering the algorithms into the various callback functions in the Matlab m-file. The steps of creating the B-RADeS GUI in Matlab are shown in the flowchart in Fig. 1 .  . This multi-level process design and simulation tool can be used to find the most efficient design of multi-stage anaerobic digester networks to achieve a defined economic and process objective. B-RADeS has several attributes that make it useful for a scientifically and economically objective process design and analysis platform for use by engineers to do their calculations during design and operation of multi-stage digesters. The main features of B-RADeS are as follows:

B-RADeS user interface
➢ B-RADeS is based on peer-reviewed models that describe growth kinetics of anaerobic digestion microorganisms. It includes ten simple biokinetic models derived based on biogas yield analogy. ➢ It does not rely on published kinetic coefficients, but it includes a section where the user determines kinetic coefficients required for digester synthesis from own experiments. Upon input of experimental data, B-RADeS automatically scans through the 10 models and ranks them in order of best fit using both quantitative and qualitative techniques. ➢ Data requirements are simple: Only experimental measurements of biogas yield are required to determine kinetic coefficients, construct attainable regions a well as synthesize digester networks. ➢ It takes into account the country-specific macroeconomic parameters (interest rate, electricity feed in tariff rate and annual working days) into the design process, which is a key motivation for investors. ➢ It is based on a systematic methodological framework for the design of multistage digester networks using the global opti- mization technique of attainable regions. The main advantage of this approach over the use of superstructure optimization is that it enables knowledge of all possible states for all possible digester configurations (even those that have not yet been devised) to be first obtained, considering mixing and biodegradation as the only fundamental process occurring in the digester.
At the level of the main interface, users can get a description of the different biokinetic models, examine the underlying assumptions and approximations of the models, fundamental concepts required to interpret AR boundary in terms of digester structures, select the level of activity (Model fitting and analysis or digerter synthesis and analysis), and export simulations results for documentation. The software interface thus allows easy interactive modeling, design and simulation of multistage anaerobic digesters taking into consideration process kinetics and economic parameters.

Digester design and analysis with B-RADeS
The following four steps are required to perform complete analysis of an anaerobic treatment process using B-RADeS: Step 1: Determine biogas yield kinetic model that best describes the organic substrate of interest. This requires data from anaerobic treatability studies using the substrate of interest. Upon input of experimental data, the software performs an automated fitting for all the ten models and ranks them in order of best fit using both numerical and graphical approaches. The numerical approach resides in the computation of a parameter, α ( Eq. (48) ), which takes into account four statistical coefficients for its computation. These coefficient include: the coefficient of determination ( R 2 ), adjusted coefficient of determination ( R 2 Adj ), root mean square error ( RMSE ) and the reduced chi-square ( χ 2 ) were major validation criteria for model selection. For good quality fit, R 2 , and R 2 Adj values should be high while RMSE and χ 2 should be low.
Models with the lowest value of α are considered more appropriate to describe a given data set if they share similar correlation coefficient. The graphical approach is based on examining the confidence contours which describe the correlation between the model parameter. The following five conditions are necessary for interpreting joint confidence regions (a) If the region, given by an ellipse is aligned with the any of the coordinate axis (vertically or horizontally), then no correlation exist between the parameters that constitute the region. (b) The parameter that lies on the coordinate axis with the greatest shadow corresponds to the parameter with the greatest variation. (c) By definition, the elliptical region is centered at the least square estimate of the model parameters. (d) If the region is a long narrow rotated ellipses, it indicates there exist significant correlation between parameter estimates. (e) If values of zero for one or more of the parameter estimates lie in the region, these parameters are plausibly zero and the corresponding terms are not significant in the model.
Models, which show less correlation between the estimated parameters are more reliable. The user can also manually test the fitting of a particular model of interest without necessarily going through the automated fitting procedure ( Fig. 3 ).
Step 2: Specify economic objective to be attained as well as county-specific macroeconomic parameters governing operation of the anaerobic digester system. The economic objective is specified in terms of the number of years required to recover investment following construction of an anaerobic digester for biogas production and electricity generation. The macroeconomic parameters are the interest rate, feed-in tariff for electricity generation from biogas as well as annual working days. B-RADeS the passes the estimated kinetic constants (for the best fitted model) to the design functionality, which together with the specified economic parameters is able Fig. 3. B-RADeS interface: Here users can input experimental data, and estimate the kinetic constants for a given digestion process. Fig. 4. B-RADeS interface: Here users set economic feasibility targets and country specific macroeconomic parameters, for B-RADeS to perform AR analysis for the user to interpret into digester structures.
to construct candidate attainable regions, and overlay the defined economic objective (payback period) onto the boundary of the attainable region in order to determine the point of intersection (see Fig. 4 ).
Step 3: Interpret the attainable regions and particularly the intersection point (which represents the optimal operating point) in terms of optimal digester structure. The interpretation of the AR boundary is based on three key fundamental results of two-dimensional AR used in everyday practice ( Ming et al., 2016 ). These include: ➢ The AR is composed of reaction and mixing surfaces only. Reaction surfaces are always convex. ➢ Points that form convex sections of the AR boundary arise from effluent concentrations specifically from PFR trajectories. ➢ Points on the AR boundary that initiate these convex PFR trajectories (from point 2 above) arise from specialized CSTRs for two-dimensional constructions.
These guidelines are provided in the Attainable region section on the main menu of B-RADeSS.

Biodigester design case studies with B-RADeS
Multi-stage anaerobic digestion in which multiple digesters are operated in a network are designed to optimize each step of the anaerobic digestion process are potentially applicable for all wastewater treatment plants ( EPA, 2006 ). Therefore, although many anaerobic wastewater treatment plants have traditionally performed anaerobic digestion processes as single stage, the use of multistate network digesters would allow these facilities to optimize the various stages of the anaerobic digestion process to meet their need. In fact, multistage digesters provide a great potential for a more efficient and flexible biogas systems that can better integrate into the bioeconomy and help harvest the energetic potential of organic waste while contributing to sustainable nutrient recycling ( Cumiskey, 2005 ;Theuerl et al., 2019 ). We illustrate the capabilities of B-RADeS considering three anaerobic wastewater treatment case studies as presented in the following section.
Case study 1: Anaerobic digestion of swine wastewater The objective here was to design a biodigester for treatment of swine wastewater that will yield a return on investment within 6 months (payback period) following start-up. Experimental data for anaerobic digestion of swine wastewater was obtained from Yang et al. (2016) . Table 4 presents fitting characteristics and kinetic coefficients for all ten biogas yield models present in B-RADeS. The models are ordered using the automating fitting functionality in B-RADeS (see step 2 in Section 4.1 ). Amongst the Ten models fitted, the Moser & Bergter, Moser and Tessier based biogas yield models were the top three based on the numerical value of the α-parameter.
However, the selection of most accurate model for description of swine wastewater required consideration of the graphical approach of confidence contours as well. Fig. 5 presents model fits and confidence contours for the first three models.
By examining the confidence contours of the top three models presented in Fig. 5 , we see that all the models have long rotated ellipses but that of the Moser-based model is narrowest indicating that there exist significant correlation between its parameter estimates. No major significant difference can be observed between that of the Moser & Bergter and the Tessier but since the former was better in terms of the numerical fitting criteria (alpha parameter), it therefore more accurately describe the experimental data for swine wastewater. Given the kinetic model of the system, the next step was to use the design analysis functionality of B-RADeS to perform an attainable region analysis of the system. Fig. 6 presents the candidate two-dimensional attainable region on which different payback period has been overlaid to indicate optimal operating points (points of intersection between the AR boundary and the objective function).
The objective function formulated is particularly important as it incorporates aspects of both total digester volume (residence time) and biogas yield. Three very interesting observations can be made from Fig. 6 . (1) For a given biogas yield, higher payback periods are achievable for larger digester volumes (higher residence times). This is because for a given biogas yield, higher digester volumes will require more investment cost, hence resulting in a longer time to break even in investment.
(2) For a given residence time on the AR boundary, shorter payback periods are achievable for higher yields of biogas. The results are highly consistent with practice because if we maintain the volume of the biogas plant constant, then we only require high yields in order to break even in investments for relatively shorter duration. (3) The range of payback periods considered intersect the AR at many points in the region, indicating that there are multiple operating points (multiple optima) for this system. However, since the objective of the design is to find a digester configuration with a payback period of 0.5 year, the reader can observe Fig. 7 , which shows how the payback period of 0.5 year has been independently overlaid onto the boundary of the AR. The left plot of Fig. 7 presents the PFR trajectory and the CSTR locus, referred to as the base trajectories.
From Fig. 6 , it is also important for readers to note that even though there are multiple intersection points of the objectives, the actual operating point to be chosen will depend on the investor's amount of capital. Points corresponding to smaller digester volume or smaller residence times (points associated with the lower part of the AR) require smaller capital investment while points corresponding to larger digester volumes require huge capital investment. Another very interesting remark in this example is that as the payback period increases, the influence of running cost (digester volume) on the payback period decreases, seen by the close proximity of the 2-, 3-and 4-year payback periods are to each other. This behavior has interesting interpretations on the investment strategy as it implies that it is more favorable to construct a larger digester, with larger operating expenses, with the intention of producing a higher biogas quantity of biogas. Hence, even if the required digester system is more complex, the plant is profitable in shorter a period of time.
Furthermore, notice from Fig. 6 that payback periods of less than 0.3 years are not achievable irrespective of the digester structure employed. Contour lines for payback periods less than 0.3 years turn to approach the horizontal axis and do not intersect the AR boundary at all. However, the attainable region is unique for a given kinetics and feed point, and if any of these change, the limits of achievability by the system may also change and hence payback periods less than 0.3 years could become attainable.   6. Candidate attainable region showing contours of different payback periods: Feed-in tariff: 17.5, annual working days: 300, discount rate: 10%, digestion time: 10days, organic load: 0.5 gVS/L, experimental methane yield: 0.5 mL/gVS. The dark green part represents a series of PFR trajectories run from each point on a CSTR locus. They are actually dark line plots but since they are many, it gives an appearance of dark green. Notice on the legend that the CSTR + PFRs appear as a dark line plot.
It is also interesting for the readers to note that the particular choice of payback period might also influence the optimal reactor structure necessary to achieve it. To attain a payback period of 0.5 year, larger capital investments will require a CSTR followed by a PFR as the optimal digester structure (corresponding to intersection point at the lower part of the AR) while smaller capital investments will require a PFR as the optimal reactor structure (corresponding to intersection point at the upper part of the AR) Case study 2: Anaerobic digestion of palm oil mill effluent (POME) Here, we consider the design and optimization of a multistage anaerobic digester for the treatment of palm oil mill wastewater in which the design objective is to attain a payback period of 1 year. We demonstrate the use of B-RADeS by considering experimental data from Faisal and Unno (2001) . Table 5 presents fitting characteristics and kinetic coefficients for all ten biogas yield models present in B-RADeS while Fig. 8 presents model fits and confidence contours for the first three models. Considering the numerical and graphical approaches for model selection (as clearly explained in case 1) the Moser based biogas yield model was selected to describe the kinetics of the process. Unlike case study 1, we notice that information for the Ierusalimsky and the Aiba based models is not included. This is because B-RADeS is programmed such that during automatic fitting of all 10 models, models that show any error during the fitting process are assigned a very large alpha value. The user then gets a displayed message stating such models and indicating that they should be deleted from the list. Hence the Ierusalimsky and the Aiba based models were not considered in the fitting experiment for POME. Fig. 9 (right plot) presents the candidate two-dimensional attainable region on which the 1 year payback period has been overlaid to indicate optimal operating points (points of intersection between the AR boundary and the objective function). The left plot of Fig. 9 presents the PFR trajectory and the CSTR locus.
Unlike case 1, the objective function intersects the lower part of the AR boundary slightly close to the feed point and rather approaches the unbounded section of the AR in the upper part of the curve. It is worth noting that the AR will always be unbounded at the residence time axis owing to the fact that states that are achieved at a given residence time will always be achievable for all later residence times ( Ming et al., 2016 ). Considering the intersection point at the lower part of the AR boundary, an anaerobic PFR is required as an optimal reactor structure to achieve a payback period of 1 year.
Case study 3: Anaerobic digestion of pharmaceutical wastewater Finally, we demonstrate the usability of B-RADeS for synthesis of a digester structure for treatment of pharmaceutical wastewater in which the objective is to attain a payback period of 2 years. The experimental data has been obtain from the work of Pandian et al. (2011) . Table 6 presents fitting characteristics and kinetic coefficients for all ten biogas yield models while Fig. 10 presents model fits and confidence contours for the first three models. Considering the numerical and graphical approaches for   9. PFR trajectory and CSTR locus (left figure) and candidate attainable region with overlaid objective function (1-year payback period). Feed-in tariff: 17.5, annual working days: 300, discount rate: 10%, digestion time: 10days, organic load: 1.0 gVS/L, experimental methane yield: 0.5 mL/gVS. model selection the Tessier based biogas yield model was selected to describe the treatment kinetics of pharmaceutical wastewater. Fig. 11 (right plot) presents the candidate two-dimensional attainable region on which the 2 year payback period has been overlaid to indicate optimal operating points (points of intersection between the AR boundary and the objective function). The left plot of Fig. 11 presents the PFR trajectory and the CSTR locus.
Unlike cases 1 and 2, the objective function intersects the lower part of the AR boundary at the feed point (0, 0), which is not feasible to operate a system at this point. However, the objective function passes through other points within the AR, any of which could be selected to operate the system. Consider a line A-B drawn such that it cuts the residence time axis as indicated on Fig. 12 .
The intersection point (point C) of line AB and the objective function is selected to be the optimal operating point and the digester structure corresponding to this point is the optimal digester design. Since this point lies on line AB, it is attainable my mixing digester effluents from points A and B, Eq. (49) (see illustration in Section 2.2.1 ). Point B lies on the PFR trajectory and is therefore attainable by running an anaerobic PFR for 3 days (point where line AB intersects the residence time axis). Point A lies on the residence time axis (corresponding to biogas yield of zero) but since a biogas  10. Model fits and confidence contours for anaerobic digestion of pharmaceutical wastewater. Fig. 11. PFR trajectory and CSTR locus (left figure) and candidate attainable region with overlaid objective function (2-year payback period). Feed-in tariff: 17.5, annual working days: 100, discount rate: 10%, digestion time: 5days, organic load: 5.0 gVS/L, experimental methane yield: 0.5 mL/gVS.  yield of zero is achieved at a lower residence time (feed point) it is also achievable at any later residence time on the residence time axis. The optimal digester structure required to achieve a payback period of 2 year therefore consist of a PFR and a bypass valve from the feed point. Table 7 provides a summary of the payback periods, intersection points as well as the required digester volumes for the three design case studies. As earlier mentioned, there a several intersection points between the objective function and the AR and the points selected in Table 7 are for illustration. In practice, the actual operating point selected by the designer will depend on other factors such as cost and space constraints. This is because different points will correspond to different digester structures some of which have different space and/or cost requirements. Fig. 13 present the optimal digester structures corresponding to the selected points of operation (see Table 7 ) for the three case studies of anaerobic digestion.
The article is of high relevance to designers of biogas digesters as it is first of its kind demonstrating the usefulness of biogas yield measurements for design and optimization of biodigester structures. Biodigester structures, which involve a staged operation of either multiple digesters or a single digester with by-pass or recycle streams has gained increasing importance due to their ability to optimize every step in the anaerobic treatment process. The authors of this study have presented a systematic model-based methodology for synthesis of biodigester structures requiring simple data requirements. The framework is based on the global optimization technique of attainable regions. The main advantage of this approach over other approaches is that it enables knowledge of all possible states for all possible digester structures (even those that have not yet been devised) to be first obtained, considering mixing and biodegradation as the only fundamental process occurring in the digester. The main novelty of the study is that it couples, biodegradation kinetics, economic objectives (payback period) and country specific macroeconomic parameters in the design process. It is also interesting for the readers to note that the particular choice of economic feasibility objective (payback period), as well as macroeconomic parameters (interest rate or feed-in tariff rate) influence the optimal biodigester structure necessary to achieve it. Due to the minimal data requirements, the study offers great promises for widespread application to enhance design of biodigester structures since biogas yield measurements are readily available from treatability studies. Even though the model-based methodolgy has been applied to only been applied to swine wastewater, palm oil mill effluent and pharmaceutical wastewater, other types of wastewaters as well as solid wastes or even wastewater sludges offers a strong research attraction for application of the framework presented in this study. This study bridges the gap between research, development and implementation of digester networks. It is interesting to compare the results of this study with our recent publication Abunde Neba et al. (2019) using attainable regions to synthesize multistage anaerobic digesters. The study considered a four-state dynamic model of anaerobic treatment process. More states imply more need for experimental measurements making it less applicable to situations where process measurements are limited. In addition, the two-dimensional attainable regions were constructed in concentration space only, and the lack of residence time makes it impossible to size the digester structure. Furthermore, the study focused on process objectives (volumetric methane productivity and gas stabilization) for design meanwhile the current study simultaneously couples process parameters, economic objectives in the construction attainable regions in residence time space, which is a key motivation for investors and makes it possible to size the digester structures.
Considering both studies put together, the results can be applied to design and optimize (based on economic and process objectives) multistage digester structures in cases of available as well as limited experimental measurements.

Conclusion
The present study was designed to develop a theoretical framework for using simplified kinetic models based on only biogas yield to model and optimize (based on economic objectives) hydrodynamic configurations of anaerobic digesters. The study has developed two-stage kinetic models based on the biogas yield approach and formulated and economic evaluation model based on simple payback period. Furthermore, the study has shown that by using two-dimensional attainable regions in residence time space, it is possible to design and optimize hydrodynamic configurations for operating low rate anaerobic digesters considering mixing and biodegradation as the only fundamental processes occurring in the digester. Attainable region analysis is a global optimization technique which incorporates elements of geometry and mathematical optimization to synthesize optimal reactor networks to achieve a given objective. For proof-of-concept, we have considered three design case studies and applied the simple payback period as the design objective for modeling optimal digester configurations.
In this article a novel software, which can be used by biodigester design engineers to rapidly model hydrodynamic configurations using experimental measurements of biogas yield. The software package has been successfully employed to model the kinetics and design optimal digester configurations for three different substrates: swine wastewater, palm oil mill effluent and pharmaceutical wastewater. Broad functionalities of B-RADeS is able to address key problems arising in design and optimization of anaerobic digester networks including: (1) modeling of anaerobic digestion kinetics by automatically fitting 10 different biokinetic models and assessing the quality of fit using numerical and graphical approaches, and finally using the selected models to determine kinetic coefficients of the process (2), Construction of twodimensional attainable regions in residence time space, and (3) optimization of anaerobic digester structures using simple payback period as well as county-specific macroeconomic parameters such as interest rate and renewable energy feed-in tariff rate. This software allows user to animate simulation results and, thereby, present them in more comprehensible and aesthetic mode. The article therefore concludes that, in principle, with only experimental measurements of biogas yield, B-RADeS can be used to generate the attainable region of the process which can be used to propose the optimal digester configuration for the process. This is highly practicable for use in small-scale onsite systems since data requirements are simple: Only experimental measurements of biogas yield are required to complete determination of kinetic coefficients, construction of attainable regions a well as synthesis of digester networks.
Finally, the study has demonstrated that the use of digester structures as opposed to single digesters improves process economics and reduces the time required to break even in investment. This result can be considered as a fundamental framework for design of digester networks using attainable regions when only biogas yield measurements are available. As recommendation for further studies, it would be interesting to apply the Infinite Di-mEnsionAl State-space (IDEAS) to obtain a general mathematical formulation for the construction of a true NRT-C-AR and compare the optimization performance with what has been obtained in the current study.

Funding
This work was supported by EnPe-NORAD under the project Upgrading Education and Research Capacity in Renewable Energy Technologies (UPERC-RET).

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.