Modeling and Dynamic Optimization of Semi-Batch Acetone–Butanol–Ethanol Fermentation with In-situ Pervaporation Membrane Separations

Abstract To alleviate the severe impact of the substrate and product inhibitions leads to the design of the semi-batch bio-reactors, which aims at finding the optimal control of substrate feed rate through time. The consequence of this problem is challenging, which highlights the numerical difficulties arising from the non-linear, constrained and often discontinuous nature of these dynamic systems. In this work, following the control vector parameterization, a profile generator algorithm is utilized to produce functions in an exponential form that can represent path of the control space comprising elements of both linearity and curvature. A novel method of stochastic dynamic optimization using Simulated Annealing incorporating direct search through a specialized model predictive controller is used to tackle final volume point state constraints. A case study of glucose to acetone–butanol–ethanol production, in terms of structured kinetics model, with in situ pervaporation will be illustrated. Prior to optimization, we proposed revised models based on model identification for the published experimental data in the literature of batch bio-reactor system. Modeling of the pervaporation membrane separations for the system have also been presented.


Introduction
Dynamic optimization (often called open loop optimal control) of semi-continuous bioprocesses is mainly aimed at finding the optimal control of process conditions (such us feed rate, temperature, pH, etc.). To identify the process, understanding insights of modeling general bioprocesses is essential. The impact of the substrate and product inhibitions, as well as the degradation of bio-substance due to metabolism inside the living cells, results in highly nonlinear dynamics with constraints present on both the control and state variables coupled with the often discontinuous nature. The consequence of this problem is challenging in terms of two issues: the construction of a profile and the development of an algorithm for dynamic optimization for the constrained model [1,2].
Clostridium acetobutylicum is a commercially applicable bacterium to produce acetone, butanol and ethanol for a century so that a fermentation process using species of Clostridium is called acetone-butanol-ethanol (ABE) fermentation [3,4]. ABE fermentation process was firstly CONTACT ta-chen Lin tc-lin@thu.edu.tw solventogenesis is in the stationary phase of cell growth; organic acids are reassimilated to form acetone, butanol, and ethanol. Thus, ABE cultivation possess complicated metabolic function resulting from substrate and product inhibitions [8]; experimentally, butanol at the level of 13.0 g/L, whereas butyric acid at a concentration of 8.7 g/L were shown completely inhibitory to the growth of cells, which lead to low productivity and yield of solvents. On the other hand, it has been shown that the specific butanol production rate could increase from 0.10 g-butanol/gcells/h with no feeding of butyric acid to 0.42 g-butanol/ g-cells/h with 5.0 g/L butyric acid (with substrate); the maximum butanol production was 16 g/L and the residual glucose concentration in broth was very low at the feeding ratio of butyric acid to glucose (B/G ratio) of 1.4. A diagrammatic illustration of the metabolic pathways of Clostridium acetobutylicum ATCC824T by Jones and Woods is capsuled in Figure 1 [6,8]. The reaction balance of target metabolites based on the structure kinetic models can be formulated as follows: patented by the chemist Chaim (Charles) Weizmann in 1916 [5] and was the primary process used to make acetone during the World War I. Thereafter, it quickly expended over the world during the early twentieth century and faded out later due to its high cost compared to petrochemicals [6]. The ABE fermentation processes using clostridia begin to regain popularity as the growing realization of bio-based fuel and chemicals on the environmental protection. Considerable researches on metabolism in species of clostridium have been carried out to develop a model describing the metabolic behaviors [7][8][9]. C. acetobutylicum is considered to have two typical phases of product formation: acidogenesis and solventogenesis [8]. ABE cultivation comprises substrate inhibition by glucose and product inhibition by acetate and butyrate in the acidogenesis phase and butanol in the solventogenesis phase [8]. For a better comprehension of C. acetobutylicum metabolism, the genome-scale model was utilized to analyze the acidogenic and solventogenic phases [8,9]. Recently, Lee et al. proposed to convert acetone to isopropanol by metabolic engineering and examine the effect of the secondary alcohol dehydrogenase on its conversion [3]. To alleviate the severe impact of the substrate and product inhibitions, fundamental researches have been proposed on fed-batch fermentation [10,11] and that with gas stripping [12] and pervaporation [4,13]; as well as, on batch and continuous fermentation [14]. Control and optimization strategies of fed-batch fermentation processes have also been demonstrated [15,16]. This article is organized as follows. Modeling of ABE metabolic pathway will be introduced first. Then a revised kinetic structured model will be discussed. In the meantime, modeling of pervaporation process and semi-batch ABE with in situ pervaporation will be demonstrated. Techniques on the estimation of control vector parameterization (CVP) and associated six parameters profile will be subsequently addressed. A case study on the incorporation of the CVP and model predictive controller into the Simulated Annealing (SA) optimization framework will be given afterward, to which most of the paper is devoted. Optimization results will be presented in the final section.

Modeling of ABE Metabolic Pathway
The metabolic pathways of ABE-producing clostridia are recognized to consist of two diverse indicative phases [8]: acidogenesis and solventogenesis. Typically, acidogenesis occurs when a cell grows exponentially, in which acetic acid (or acetate) and butyric acid (or butyrate) are produced with ATP formation. The acetate and butyrate are toxic to the cell at high concentrations. Subsequently, Since many metabolic reactions in ABE production occur in the presence of ATP or NADH, these reactions may terminate in case of insufficient energy, i.e. after the exhaustion of glucose. In regard to the cessation of metabolic reactions after glucose exhaustion, an on-off mechanism was introduced into the rate equations of specific metabolic reactions (r1 to r19) in Figure 1. Introduction of on-off mechanism can be referred to Model III in the reference by Shinto et al., whose kinetic parameters were estimated to capture the experimental time-course data, where the r1 and r19 are shown in the paper as follows [8]: where F s is value 1 or 0; its value is assumed to be 1 if glucose concentration is over 1.00 mM and 0 if it is under 1.00 mM. Here, we further clarify the r1 and r19 based on the simulation results shown in the paper: [Byturate] ⋅ F s notes: enzymes are indicated in bold and abbreviated as follows: Pta, phosphotransacetylase; aK, acetate kinase; coat, coa transferase; PtB, phosphotransbutyrylase; BK, butyrate kinase; BaDH, butyraldehyde dehydrogenase; BDH, butanol dehydrogenase [6,8].
diagram of incorporating the model identification into the SA optimization framework. To solve initial value problems for ODEs, the algorithms in the ode45 solvers is used for 'nonstiff ' problem type in the medium order of accuracy. Here, the possible range of decision parameters d 1 -d 4 and f 1 -f 4 are all [0, 600]; however, if the number of the parameters in each inhibition terms in Equation (22) being further reduced to three parameters, then the upper bounds of the possible range can be estimated to be smaller. Optimal kinetic parameters are illustrated in Table 1 by fitting Equation (22) to experimental data of 70.6 mM of initial glucose with the value of SSE reduced from 4.4708 (original model as in Equation (21)) to 2.2955 (new model as in Equation (22)). As shown in Figures 3 and 4 (with the same kinetic parameters in Table 1 but 122 mM of initial concentration of glucose), the simulation results showed much more qualitative consistency with the experimental time-course data.

Modeling of Pervaporation Process
The mathematical model of Wenchang and Sikdar was utilized to simulate the pervaporation process [18]. The overall resistance is constructed by the three phases, namely, 'the bulk, the membrane and the vapor' as illustrated in Figure 5. It is assumed that the mechanism describing the mass transfer of the solute consists of five steps [13]. This includes: (i) diffusion from the bulk feed stream to the feed-membrane interface, (ii) absorption into the membrane, (iii) diffusion through the membrane to the downstream, (iv) desorbtion into vapor at the permeate side, (v) diffusion from vapor-membrane interface to the bulk of the vapor phase.
Thus, the permeation flux through the boundary layer on the feed side of the membrane is formulated as: where k L is the liquid boundary layer mass transfer coefficient; ρ L is the total molar volume concentration in liquid at the feed side; X ib and X * i are the mole fractions of the component (i) at the bulk feed stream and the membrane interface, respectively. In this pervaporation model, superscript '*' and subscript 'b' denote interface and bulk, respectively. If P * i is the partial pressure at the membrane interface at the permeate side and Q m i is the vapor permeability of the component (i) through the membrane (mol m/m 2 kPa s), the flux through the membrane may be expressed as: The rate equation −r13 in Equation (7) indicating the death reaction of the cell is formulated in the reference [8] as follow: We further propose to modify the rate equation r13 as in Equation (22), in addition to original model whose values of all the kinetic parameters of 70.6 mM of initial glucose was used [8]: It should be mentioned that when glucose concentration is high compared to constant f 2 and f 4 , the glucose inhibition term will be reduced to ; similarly, butanol concentration is high compared to constant d 2 and d 4 , the butanol inhibition term . Note also that the four parameters d 1 -d 4 or f 1 -f 4 can be further reduced to three parameters in case of all of the four parameters divided by selecting one parameter for d i The optimization is implemented using SA Algorithm [17], which is a derivative-free stochastic method capable of approaching the global optimum. This assumption of modified model in Equation (22) is valid to minimize the sum of square error (SSE), as the objective function (G) formulated generally in (23), between experimental data (X exp hik ) and simulation results (X pred hik ).
where where X pred is the vector of predict state variables of the ordinary differential equations (ODEs) as expressed in Equation (24) and p is the vector of independent or decision parameters. In addition, nh, ni, and nk are the number of the experiments, the state variables and the sample time points of the experimental data, respectively. It is solved by the 'simulannealbnd' solver in Matlab R2010a software package, which aims to find unconstrained or bound-constrained minimum of function of several variables using SA algorithm. Figure 2 demonstrates the phase boundary layer, where k V is the mass transfer coefficient in the vapor boundary layer and P ib is the partial vapor pressure at the permeate side, may be defined as: where, for the dilute organic aqueous solution, Here, l is the membrane thickness, γ i is the activity coefficient, ∞ i is the infinite dilution activity coefficient, and P o i is the saturated vapor pressure. The flux through the permeate vapor  Reaction  pervaporation cell. Therefore, the change of volume inside the tank (V) with time can be formulated as follows: where A m is the membrane area (cm 2 ). The change of the component (i) molar mass with time will be: By substituting Equation (29) in (30), the change of the tank concentration with time can be described by: The experimental results, as shown in Figure 6 where the acetate and butyrate concentration are not measured, for poly(ether-block-amide) (PEBA) membrane in the reference [20,21] are used to verify the developed model as Equations (28) and (31). The change of end-products concentrations in the feed tank during the pervaporation process was predicted using the developed model. The pervaporation membranes PEBA were synthesized with a total active area of 0.08 m 2 and an average thickness of 50 μm. To investigate the pervaporation performance, a model solution consisting of 4 g/L of acetone, 8 g/L of butanol, 1.33 g/L of ethanol, 3.5 g/L of acetic acid and 3.5 g/L of butyric acid was prepared to imitate the supernatant of cultivation broth. To decide the optimal parameters (K i and β i ) in Equation (28), SA Algorithm is implemented to minimize the SSE between experimental data and simulation results. Here, the objective function is the value of SSE and the possible ranges of all the At steady state conditions, the molar flux has the same value at any location in the membrane. This may lead to a simple flux expression of Equation (26), by combining Equations (23)- (25), in which it is related directly to the downstream pressure.
where Q i is the overall permeability [18] and K OV, i is the overall mass transfer coefficient across the membrane: To further simplify the molar flux for the modeling of the unknown parameters, Equation (26) becomes as follows: where K i is overall mass transfer coefficient (m/s), which is equal to K OV,i α i ; here, K i is revised from El-Zanati et al. [13]; C ib and ρ L are concentration of the component i and total component, respectively, at the bulk feed stream. Therefore, L = ∑ C ib such that C ib = ρ L X ib . The concentration C ib is calculated as the mean value between the start and the end feed concentration, which is also written as C bm [11,13] . In case downstream pressure is very small, 'P ib~0 ' , (28) can be simplified as J i = K i C bm . If the change of the concentration with time is linear, the arithmetic mean concentration is used, while logarithmic mean concentration can be used if it is nonlinear [19]. The change of mass of the contents of the feed tank is actually equivalent to the amount permeated by the total molar flux (J tot = ∑ J i ) or the total flux (g/cm 2 s) via the (28) Figure 4. experimental time-course data and simulation results of 122 mm of initial concentration of glucose.
As indicated in Figure 4, the five end products of the metabolic pathways include: butanol, ethanol, lactate, acetate, butyrate, and acetone, which can be removed by membrane of PEBA. Therefore, batch models in Equations (8), (9), (13), (14), and (16) can be modified by the generation described in Equation (35). In addition, the only state constraint considered is for the final volume as follows: Here, V T is the volume of the tank, which is equal to 300 mL.

The Six-parameter Profile Generator Algorithm
Dynamic optimization of the semi-batch process is intended to find the optimal feeding profiles through time.
In order to apply the SA algorithm to dynamic optimization of semi-continuous bioprocesses successfully, the construction of a control profile is required that mirrors subtle nature of biological kinetics. An appropriate methodology must allow discontinuous operations and nonlinear curves. In view of the complexity of CVP [22,23], an efficient mechanism of generating various control profiles that approximates the basis of the optimal control space is vitally important. The profile generator algorithm proposed by Choong [24] has shown to capture the comprehensive features of different curves comprising concavity, convexity and linearity by using six parameters. Basically, the mathematical formulation of six-parameter profile is optimization parameters are [0, 5]. Table 2 illustrates the optimal decision parameters of the above pervaporation model for 10 executions of the program codes. Figure 6 depicts the experimental data and predicted concentration in feed tank by the developed model with optimized parameters.

Modeling of Semi-Batch ABE with In-situ Pervaporation
To develop models for semi-batch ABE fermentation with in situ pervaporation, the volumetric feeding rate (F) of glucose is introduced into the batch process formulations. Thus, the change of glucose and biomass concentration with time in Equations (1) and (7) become as follows: For the end products of metabolic pathways in batch mode, the change of end products concentration with time can be generalized as: where r p is the lump of bio-kinetics rate equations contributing to products formation. As mentioned in the pervaporation by membrane of PEBA, the permeate flux results in the loss of product concentrations as illustrated in Equation (31). Therefore, we reconstruct the general model for the change of products concentration with time as the following: where the initial conditions of both integrals is assumed to be zero and the solution of the integration constant in Equation (39) can be derived. If the rate of the feed is described by a six-parameter profile, the constraints imposed on the total volume charged into the bioreactor will be related to the integration of the six-parameter profile. This may incur feasibility problems where the constraints imposed on the control and state variables cannot be applied explicitly, and the solutions can only be directly controlled by implicitly the estimated feasible bounds of decision parameters.

Dynamic Optimization using SA
Dynamic optimization using SA is in fact the combination of simulation of processes models coupled with optimization framework in two separate domains [1]. The framework of SA involves the transition of objective value or performance index generated by randomly perturbing states or control variables based on processes of the Markov chain. As shown in Figure 8, the only interface between optimization and simulation processes is the evaluation of objective function. Here, the decision parameter for perturbation including the following: ; a 1 and a 2 ≥ 1, t total ∈ [0, 200]; T on and T off are the time when the in situ pervaporation starts and terminates, respectively. Here, the practice of turning on and off of the pervaporation is by adding the relationship of Am = 0 if t is outside the period of [T on , T off ]. During the course of the optimization, the various profiles in an individual unit are governed by the above decision parameters; in each loop, a random number generator is used to produce a real number between the between 0 and 1. This real number is then compared against various perturbation probabilities. The choice of the control variable for perturbation is random and dictated by its perturbation probability. If the perturbation probabilities for the control variables are equal, then the total perturbations implemented at the end of the optimization are the same for each control variable. The perturbations of the chosen control variable are used to produce a real number between the user-specified lower and upper bounds.
[X 0 , X F , X inter , tb 1 , a 1 , a 2 , t total , T on , T off , Am] composed of two distinct types of exponential function in the form of: Type 1: where the function X(t) is defined in the interval of operating time t ∊ [0, t total ] and the subscript 0 and F denote the initial (at t = 0) and final values (at t = t total ), respectively. In addition, the power of the exponential function a 1 and a 2 are supposed to be positive numbers and greater than one. It is obvious that the difference between these two types of exponential function lies in its base variables, such that the power a 1 and a 2 represent convexity and concavity of the curves relatively. In addition to the above four parameters, the profiles can be controlled by connection of the two curves: the conjunction value (X inter ) and the scaling time location (tb 1 ). The value of 'X inter ' indicates the conjunction point of the two functions. The absolute value of tb 1 denotes the division of subinterval where the two curves join, while the sign of the value of tb 1 implies the arrangement of the type sequence. The positive intermediate peak location refers that the Type 1 function being used in the first sub-interval of 0 < t < abs(tb 1 ) and the Type 2 in interval of abs(tb 1 ) <t < 1. The negative intermediate peak location generates the opposite profiles [25]. Figure 7 illustrates the shape of profiles of the (I) Type1 + Type2 and (II) Type2 + Type1 corresponding to positive and negative intermediate point location.

The Integration of the Six-parameter Profile
In order to predict the liquid volume accumulated in the bio-reactors, the feeding flowrate constructed by the six-parameter profile is integrated mathematically. The analytical formulation of the integration of the six-parameter profile comprises the integral of Equations (37) and (38), respectively, as follows [24]: Type 1: control the feasible range for decision parameters within the range between two extreme bounds. However, the extreme bounds of the decision parameters that satisfy the optimization constraints relating to state variables may not be known explicitly. In addition, it is often the case that the global optimum exists at the boundaries of the constraints. Therefore, an alternative technique of MPC can be employed to solve the kind of indirect bounds for decision parameters subject to dependent infeasibility constraints.
The MPC is incorporated in order to predict the next feasible control corresponding to controls and state trajectories constraints, respectively. As can be observed in Figure 9, the MPC controller is employed in terms of two coordinates: the iteration of the new profile generation in the optimization framework (vertical iteration axis) and time domain in simulation processes (horizon time axis). At every entry of a new iteration point (i k ) in the SA optimization algorithm, one of the decision variables is selected according to the probability. The chosen decision variable is then randomly perturbed between the predefined upper bound and the lower bound. Thus,

Incorporation of Direct Search via Model Predictive Controller
The model predictive controller (MPC) is incorporated in order to predict the next feasible control corresponding to controls and state trajectories constraints, respectively. In addition to trial and error, several controlled search strategies can be exerted on the profile during the course of the optimization to tackle the feasibility constraints on both controls and state variables. One class of the major controlled search methodologies is the use of search region reduction and direct search methodology. In the former, the search region reduction can be implemented by reducing the sampling distribution according to the distance moved away from current value. The objective of the direct search methodology is to predict and   Table 4. Optimal decision parameters for the semi-batch aBe with in situ pervaporation (S 0 = 70.6 mM).

Optimal results
Decision parameters     Table 5. Optimal decision parameters for the semi-batch aBe with in situ pervaporation (S 0 = 150 mM).

Decision parameters
Optimal results

Numerical Study
The problem formulated in the previous section results in a dynamic optimization using SA problem. It is also solved by the 'simulannealbnd' solver in Matlab R2010a software package, to solve initial value problems for ODEs, the ode45 solvers is used; however, we do notice that when the choice of the different ODE solvers can tackle 'stiff ' instant, t k , which starts from the initial point, and moves into next sampling instant in the direction of the end of operating time. The model predictive controller works in the same manner as the integrator but the controller is designed to minimize error due to the bias between processes outputs (i.e. state trajectories) and the desired set points.  glucose concentration of the feeding substrate is set to be 500 g/L and thus, with density of glucose 0.18 g/mmole, the total volume of the feed can be calculated. Table 3 lists the initial concentrations of each metabolite in the semibatch ABE fermenter with in situ pervaporation. Tables  4 and 5 show the optimal decision parameters results for the semi-batch ABE fermenter with in situ pervaporation problems when the control profiles appear to be discontinuous. Here, the objective function is the equivalent total moles of butanol production, which includes the butanol remained in the bioreactor and the butanol separated by the pervaporation membrane. The design of the fed-batch cases are structured to compare to the batch experiments using the equivalent quantity of substrate. Here, the The dynamic optimization results indicate that, with the same quantity of glucose used in batch processes, the optimized substrate feeding profile coupled with in situ pervaporation could improve the total butanol production to be equal to 22 corresponding to initial glucose concentration (S 0 ) used in 70.6 mM and 150 mM batch bioreactors, respectively. The best optimized fed-batch with pervaporation is illustrated Figures 10 and 11, which are equivalent to glucose used in 70.6 and 150 mM batch bioreactor, respectively. The latter case is compared to 122 mM batch experimental results referred to reference [8], which is equal to 28.5 mmole. It should be mentioned that Figure 10 illustrates the best performance index of total butanol 22.5897 (mmole) in Table 4, where the total reaction time (T total ) is 192.0147 (h) and the time when the in situ pervaporation starts (T on ) and terminates (T off ) are 178.5962 (h) and 182.2896 (h) with Am being 0.0238 (m 2 ). As shown in Figure 10, the flow rate terminates due to the full reactor volume at about 120 (h), and afterward the butanol concentration reach its highest concentration and therefore the pervaporation starts. In addition, Figure 11 represents the best performance index of total butanol 41.8452 (mmole) in Table  5, where the total reaction time (T total ) is 179.6663 (h) and the time when the in situ pervaporation starts (T on ) and terminates (T off ) are 132.4463 (h) and 135.4117 (h) with Am being 0.0045 (m 2 ). As shown in Figure 11, the flow rate terminates due to the full reactor volume at about 136 (h) , and however the butanol concentration reach high concentration so as to activate pervaporation before terminal feeding. In the above-mentioned two cases, the duration of the pervaporation time are very short compared to total reaction periods, where perhaps the low flowrate and long period of reaction have compensated the inhibitions of products.

Conclusion
In this study, we propose the modified rate equation indicating the death reaction of the cell by introducing glucose and butanol inhibition terms. Secondly, the model identification of the modified kinetic and pervaporation model parameters using the SA algorithm. Finally, dynamic optimization incorporating into the SA using direct search by the specified model predictive controller is demonstrated. A case study of glucose to ABE semi-batch production with in situ pervaporation has been illustrated to improve final butanol produced comparing to batch processes. In the future, if the simulation results could further be validated by the experimental data, the intensification on specific region of process design would be possible by reconstructing the precise models to enhance the practical applications.