Closed-Loop Control of Fed-Batch Bioreactors: A Shrinking-Horizon

Ind. Eng. Chem. Res. 2004, 43, 3381-3393

3381

Closed-Loop Control of Fed-Batch Bioreactors: A Shrinking-Horizon Approach Abhishek S. Soni and Robert S. Parker* Department of Chemical and Petroleum Engineering, University of Pittsburgh, Pittsburgh, Pennsylvania 15261

The fed-batch mode of bioreactor operation is commonly used in the chemical and biopharmaceutical industries. Bioreactors typically are characterized by highly nonlinear behavior occurring on both a macroscopic reactor scale and a microscopic cellular scale. This nonlinearity, coupled with the dynamic nature of the fed-batch mode, makes for a challenging control problem. In this work, a multiscale model1 describing the growth of yeast in an aerobic fed-batch mode is employed to address both the offline and online optimization and control issues for the fedbatch fermentation problem. Nonlinear optimization techniques are used offline to generate an optimal substrate feed profile for the nominal problem that maximizes the end of the batch ethanol concentration; shrinking-horizon nonlinear quadratic dynamic matrix control is used online for closed-loop trajectory tracking and disturbance rejection. Nominal performance is admirable, resulting in final ethanol concentrations within 1% of the open-loop optimal ethanol concentration. An online re-optimization strategy is presented for disturbance compensation, which improves the final ethanol concentration in the presence of disturbances by 17% when compared to the case when no disturbance compensation strategy is employed. 1. Introduction Rising energy costs, increased global competition in terms of both price and product quality, and the desire to produce in a “green” manner have paved the way for a biological route toward manufacturing. An additional consideration is that many of the products obtained by the biological route, such as monoclonal antibodies, proteins, and other therapeutic drugs, either cannot be produced or are very difficult to obtain by conventional manufacturing methods. Most of these are classified as low-volume/high-value products, and there is a tremendous economic potential in this sector of the market. However, the trend in these processes is to operate in a so-called (fed-)batch mode to maintain sterile conditions, avoid substrate inhibition, and allow operational flexibility. This operation generally follows a set pattern akin to a recipe and is repeated from one batch to the next. These recipes are usually heuristically developed by a trial-and-error procedure requiring significant experimental effort.2 By using the available biological information in combination with tools from control theory, it is possible to develop effective control strategies that could avoid some of the empiricism associated with the operation of (fed-)batch bioreactors and, in turn, lead to reduced product variability and optimal resource utilization. Bioreactor operation is inherently complex owing to the reactions occurring at both a macroscopic reactor scale (such as substrate uptake, oxygen transport, etc.) and a microscopic cellular scale (described by the metabolic pathway of the cell).3 The equations that are used to describe these processes, and thus the system model itself, are generally nonlinear. In addition, the fed-batch mode of bioreactor operation is dynamic in nature with no single steady state during the course of the batch. Furthermore, the process operation is usually * To whom correspondence should be addressed. Tel.: (412) 624-7364. Fax: (412) 624-9639. E-mail: [email protected].

constrained by either product quality or operational requirements. All of these factors combined make for a challenging control problem. Because of its economic and industrial relevance, the fed-batch bioreactor control problem has received considerable attention. Open-loop control involves the determination of a substrate feeding policy that maximizes the product concentration in the absence of disturbances, and this problem has been solved by a number of different techniques. Chen and Hwang4 reported an on-off control scheme to obtain the optimal substrate policy, maximizing the final ethanol concentration for substrate- and product-inhibited fed-batch baker’s yeast fermentation. The input was parametrized over the batch period, and the problem was solved by transforming it to a nonlinear programming (NLP) problem. Dynamic programming has also been used to solve the dynamic optimization problem.5-7 Tholudur and Ramirez utilize a neural network model to predict the functional forms (as opposed to state information) to model the bioreactor dynamics. This neural network parameter function model is used in conjunction with iterative dynamic programming5 to solve for an optimal profile. Pushpavanam et al.8 considered a sequential quadratic programming approach in which the entire batch was divided into equally spaced intervals or stages and the feed was introduced as discrete pulses at the beginning of each interval. Banga and co-workers9 used a control vector parametrization technique and transformed the original optimal control problem (OCP) into a NLP problem. They presented an algorithm in which deterministic elements of the solution such as Hessians and gradients for the objective function were calculated using second-order sensitivities. In contrast to the approaches considered above, efforts have also been directed to obtain the control policy in a feedback or closed-loop form. Lim and co-workers10,11 were among the first to analyze the problem of feedback optimization for the singular control problem of fedbatch fermentors. They formulated the problem based

10.1021/ie030535b CCC: $27.50 © 2004 American Chemical Society Published on Web 05/21/2004

3382 Ind. Eng. Chem. Res., Vol. 43, No. 13, 2004

machinery. Closed-loop control approaches have been primarily limited to expressing the solution in a feedback form, and online control implementation has not received widespread attention. In this work, a multiscale model is used, which can capture the important phenomena of diauxic growth, i.e., sequential uptake of glucose and ethanol as substrates, in fed-batch cultures of baker’s yeast (Saccharomyces cerevisiae). The term multiscale in this work is with reference to the size scale rather than the time scale and results from including considerable biochemical detail at the cellular level without resorting to the use of another “cellular” time scale. Using the multiscale model, open-loop optimal feed policies are first determined. An online closed-loop controller is then developed in order to track the offline optimal trajectories with a special focus on strategies for disturbance compensation. Thus, the focus of this work is on the practically relevant issue of online closedloop constrained control of fed-batch bioreactors, an area that is largely unexplored.9 Figure 1. Coupling of microscopic and macroscopic models through the specific growth rate µ.

on singular control theory and obtained a nonlinear relationship between the feed rate and the system state variables in a feedback form. Palanki et al.12 analyzed the OCP for baker’s yeast fermentation from a geometric perspective and derived optimal feedback laws for the singular region of operation. Hodge and Karim13 utilized a nonlinear model predictive controller (NLMPC) for fed-batch fermentation of Zymomonas mobilis, with control and prediction horizons of unity so as to approach the setpoint in a single control move. Most of the work on fed-batch fermentation control has focused on the open-loop OCP. Furthermore, the openloop control strategy has been demonstrated using models with a simplistic description of the cell’s biochemical

2. Bioreactor Process Modeling Because bioreactor operation involves actual living cells, an adequate representation of the bioreactor operation must include descriptions of events occurring at both the reactor and cellular scales. In this regard, one of the key variables that describes microorganism growth in a bioreactor model is the specific growth rate, µ. In the multiscale model employed here, µ serves as the coupling parameter between events occurring at the macroscopic and microscopic scales; a general schematic of a multiscale model coupling is shown in Figure 1. At the macroscopic level, the system is described by mass balances for the biomass, substrate(s), reactor volume, and product, whereas the biological microscale component is a combination of the metabolic pathway and the kinetic information. The specific growth rate and prod-

Figure 2. Schematic for the distribution of the central glucose flux in S. cerevisiae.

Ind. Eng. Chem. Res., Vol. 43, No. 13, 2004 3383

uct formation rate are dependent on the rate of substrate uptake from the bulk medium surrounding the cell and on the rate of reaction in the cellular metabolic pathway. The process under consideration in this work was a fed-batch bioreactor growing S. cerevisiae on glucose. The model was described by Enfors and co-workers1 and is summarized here for completeness. The process model includes a description of the metabolism in S. cerevisiae, and this is used for the simulation of concentrations of the state variables, biomass (X), glucose (S), ethanol (E), medium volume (V), and dissolved oxygen tension (DOT). The rate equations for the metabolic fluxes are expressed as specific rates q (g/g‚h), and the notation used here is the same as that used in the original paper. A schematic for the central glucose flux distribution is shown in Figure 2. The first step in the flux pathway is the substrate uptake (qS), which is assumed to follow Monod kinetics and includes an experimentally observed lag-phase term. This flux is then distributed into the oxidative (qSox) and fermentative (qSf) pathways. The oxidative flux is, in turn, composed of fluxes for anabolism (qSoxan), maintenance (qm), and energy for growth (qSoxen). The flux of oxygen needed for sugar oxidation (qOS) is a function of the stoichiometry of respiration and qSoxen. However, qOS is limited by a maximum respiration capacity (qOmax) so that when qOS > qOmax, there exists a proportional reduction in the fluxes of qSox and qSoxen. One of the unique features of the model is that qOmax is not constant but is a function of the instantaneous DOT and ethanol concentration. The part of qS used for the fermentative pathway is obtained by a difference with the knowledge of qS and qSox. As with aerobic metabolism, qSf is split again into fluxes for anabolism (qSfan) and for energy metabolism (qSfen). The flux for ethanol production (qEp) is then obtained from the stoichiometry of the reactions from glucose to ethanol. The model also includes a description of ethanol consumption (qEc). The consumption of ethanol as a substrate by the cells occurs during aerobic cultivations when the glucose flux is low such that qOS < qOmax. This rate of ethanol consumption is, in turn, utilized for anabolism, qEan, and for energy utilization, qEen. The flux for oxygen uptake (qO) is also dependent on qEen. Finally, the specific growth rate µ is calculated as a function of qSox, qSf, and qEc; this includes growth contributions from using glucose and/or ethanol as substrates based on the system operating conditions. Using these rates, mass balance equations for the state variables are given as14

dX F ) µX - X dt V

(1)

dS F ) (S - S) - XqS dt V i

(2)

dV )F dt

(3)

dDOT ) KLa(100 - DOT) - 14000XqO dt

(4)

dE F ) X(qEp - qEc) - E dt V

(5)

The equations for biomass, substrate, and ethanol consist of a positive growth or accumulation term and a second term containing the factor F/V due to the

dilution of these components as the feed is added during the reactor operation. The concentration of substrate in the feed is Si, and it is assumed that the feed is sterile; i.e., it does not contain any biomass. KLa is the liquidside mass-transfer coefficient in the balance equation for dissolved oxygen. Changes in the model equations and in the sequence of calculations from the original paper1 are reported in the Appendix.14 2.1. Open-Loop Simulation Results. Open-loop simulation runs were carried out to characterize the system behavior as a function of changes in the substrate concentration. The state variable concentration profiles were obtained for a batch time of 20 h. The initial conditions were selected so that the final volumes for both the batch and the fed-batch reactors were 20 L. Also, both reactors were supplied the same mass of glucose over the course of the batch. For the fed-batch case, a sterile feed was selected with concentration Si ) 100 g/L (for comparison with refs 8 and 9). Model simulation results of batch and fed-batch growth on glucose are shown in Figure 3. The solid lines show the concentration profiles for batch growth with an initial substrate concentration of 100 g/L. During the initial growth phase, glucose consumption leads to biomass growth and ethanol production with a fairly high specific growth rate. As the batch continues, the specific growth rate decreases until all of the initial glucose charge is exhausted. Biomass growth then continues using ethanol as a substrate, as a result of which the concentration of ethanol decreases. The dash-dotted lines represent the concentration profiles for a fed-batch mode of operation with a specified input profile as the feed (Si ) 100 g/L). For the batch phase of the fed-batch fermentation, the specific growth follows a trend similar to that observed during the batch growth. However, once feed is introduced, the specific growth rate increases because of the glucose available from the feed. As the glucose from the feed is consumed, the specific growth rate decreases. The availability and continued consumption of glucose ensures that ethanol is not consumed as a substrate. Thus, the model is able to capture the dynamics of sequential growth on glucose and ethanol in both batch and fed-batch mode. If ethanol is intended to be the product of interest, then the consumption of ethanol is undesirable. Unlike the model considered above, models that have been previously used to describe fed-batch baker’s yeast fermentation8,11 do not include a detailed representation of the cellular level growth process, mainly the flux for ethanol consumption. As a result, they are unable to capture the diauxic growth phenomenon observed in (fed-)batch cultures of baker’s yeast. With the knowledge that the model can capture the consumption of ethanol, offline optimization can be carried out to determine the feed profile that maximizes the ethanol concentration at the end of the batch. 3. Open-Loop Optimal Control For the fed-batch bioreactor studied in this work, the goal is to maximize the final ethanol concentration. Hence, it is required to determine the input profile that would maximize the end-of-batch product concentration, and this forms the conceptual basis for the OCP. 3.1. OCP: General Formulation. Consider the system dynamics described by

x3 ) f(x,u,t)

for t g t0 where t0 is fixed

(6)


Figure 3. State and specific growth rate profiles in batch (s) and fed-batch (-‚-) growth conditions, the latter with the glucose feed profile (- -). Initial and operating conditions are given as follows. Batch: X0 ) 0.01 g/L, S0 ) 100 g/L, E0 ) 0 g/L, and V0 ) Vf ) 20 L. Fed-batch: X0 ) 0.2 g/L, S0 ) 100 g/L, E0 ) 0 g/L, V0 ) 1 L, Vf ) 20 L, and Si ) 100 g/L.

In this equation, x(t) and u(t) are the vector-valued state and input, respectively, and t0 is the initial time. Associated with the process operation is an objective function that needs to be maximized, and the general formulation is given as

J(u) ) φ[x(tf),tf] +

∫tt L(x,u,t) dt f

0

(7)

The function φ accounts for the contribution of the final state, whereas L accounts for the path dependence in the objective function with tf as the final time of operation. Additional constraint equations may need to be included, and they are usually present in the form of bounds on the state and/or manipulated variables. The detailed solution technique for the OCP can be found in ref 15, and only the main steps in the solution procedure are presented here. Initially, the necessary condition for u to be optimal is that it should maximize the Hamiltonian, which is obtained by appending the system dynamics to the objective function using a costate variable so that

H(x,u,t) ) L(x,u,t) + λTf(x,u,t)

(8)

In this equation, λ is the co-state variable and is used to incorporate the system dynamics. The original OCP is then transformed into a two-point boundary value problem because the differential equations for the state and new co-state variables have boundary conditions defined at t0 and tf, respectively. 3.2. Solution Approaches. Different approaches have been used for the solution of OCPs for fed-batch fermentors, and these can be broadly divided into two distinct categories, one based on Pontryagin’s maximum principle11,16 and the other based on nonlinear optimization techniques.8,9,17 When the substrate feed rate is used as the control variable for the fed-batch fermentation problem, the Hamiltonian is linear with respect to the control variable. Consequently, the solution of the OCP is either

bang-bang for nonzero values of the Hamiltonian or singular when the Hamiltonian is zero. In the singular case, both the Hamiltonian and its first and higher derivatives are zero so that the maximum principle does not provide any information for determining the optimal control input. As a result, part of the optimal trajectory may be singular, and switching times between the singular and nonsingular regions need to be determined. Differential equations for both the state and co-state variables then need to be solved in the different regions, and the solution has been reported to be highly sensitive to the choice of the initial conditions and the initial guess values used for the solution of the optimization problem.8 In addition, when constraints are present on the control variable and/or the system state variables, the resulting problem poses numerical and computational challenges. The applicability of the maximum principle is therefore limited to low-order problems, or to problems with fairly low complexity in the process model. As an alternative to Pontryagin’s maximum principle, the original OCP can be converted into a NLP problem. The mathematical basis for this solution technique is the discretization of the original continuous input profile in order to make the problem computationally tractable. In this method, the fed-batch interval (t0, tf) is partitioned into subintervals and the control variable is approximated by suitable functions, such as piecewise constant functions, within each interval. The resulting problem is thus reduced to a mathematical programming problem that can be solved using standard optimization algorithms. The switching times are pre-assigned whereas the heights of the piecewise constant functions are decision variables for the optimization algorithm. 3.3. Bioreactor Case Study. For optimal control of the fed-batch ethanol fermentation considered in this work, the state variable x and the system equations f are described by the multiscale model represented by eqs 1-5, so that

x ) [X S V DOT E ]T

(9)


Figure 4. Plot of state variable profiles for open-loop optimal feed delivery. Left: biomass (s), volume (-‚-), and optimal input (- -). Right: glucose (s) and ethanol (-‚-). X0 ) 0.2 g/L, S0 ) 100 g/L, V0 ) 1 L, E0 ) 0 g/L, and Si ) 100 g/L.

The feed rate F serves as the control input u, and the objective function is to maximize the ethanol concentration at the end of the batch, given as

max J[u(t),x(t)] u(t)

J ) E(tf)

(10)

In addition, it is assumed that the system is constrained so that the reactor volume at the end of the batch is restricted to 20 L, which is an end-point constraint. Furthermore, the feed rate cannot be negative, and the feed rate integrated over any interval cannot exceed the free reactor volume at the start of the interval. The final batch time was chosen as 20 h, and the entire batch was split into equal intervals of 1 h duration. The control input was then discretized into 20 piecewise constant segments with values u(k) for k ) 0, 1, 2, ...,19, where k is the time instant at which the input is initiated. This problem formulation results in the following feed rate and volume constraints:

0 e u(k) e 20 - V(k)

(11)

19

∑ u(k) ) 20 - V0

(12)

k)0

In these equations, V(k) is the volume at time k and V0 is the initial volume. The optimization problem was solved using the routine fmincon in MATLAB (The MathWorks, Natick, MA), and the results of the optimization strategy are shown in Figure 4. The optimal control policy for the fed-batch fermentation shows two distinct regions. Initially, it is seen that the fermentation proceeds in a batch mode when the value of the substrate feed rate is zero; substrate addition then follows the batch period. The batch growth continues until all of the initial substrate is completely consumed. To prevent the consumption of ethanol as a substrate, the glucose feed is sent into the reactor and is utilized for the production of ethanol in preference to biomass. This approach leads to the optimal glucose feed profile for maximizing the end-of-batch ethanol concentration. The results of this approach are qualitatively similar

to the optimal profile obtained by Pushpavanam et al.8 and by Modak and Lim;11 i.e., the optimal profile is first batch in nature, followed by a fed-batch mode. Quantitatively, the final ethanol concentrations obtained are different because they are obtained by generating optimal profiles from different model structures. The offline optimum profile provides the reference trajectory that the fed-batch operation must follow to maximize the end-of-batch ethanol concentration in the absence of disturbances. The optimal control policy represents the best performance that can be obtained from the system for the given set of initial and feed conditions. 4. Closed-Loop Control Most of the control literature for fed-batch cultures focuses on open-loop operation owing to the highly nonlinear and inherently difficult dynamic behavior.18 Optimization is carried out offline, and the reactor is fed with the determined optimal feed profile. Once the batch proceeds, there is no provision to account for disturbances that occur during the batch. In the case of run-to-run control,19 the parameter values used in the system model are updated at the end of the batch based on the progression of the previous batch. On the basis of the new values of these parameters, reoptimization is carried out and the new optimized policy is applied to the next batch. Thus, the problem of implementing the optimal policy online and ensuring that the system follows the optimal trajectory in the presence of disturbances has received limited prior consideration.9 Attention is thus focused on designing a controller to track the optimal policy and on approaches for disturbance compensation for the closedloop control problem. Conventional MPC,20-24 which employs a recedinghorizon framework, is geared more toward continuous processes.25 In contrast to continuous MPC, the control horizon in batch processes shrinks as the batch nears completion. Thus, the shrinking-horizon MPC (SHMPC) formulation is utilized in this study because it is more suited to batch processes, where the available window for control shrinks as the end of batch nears. The SHMPC approach has been utilized by Liotta et al.26 for control of the particle size in a semibatch emulsion


polymerization reactor and by Thomas and co-workers for quality control of composite laminates.27 Considering NLMPC, a nonlinear optimization problem needs to be solved at every time step. The general methodology is to convert the objective function and the constraints into a NLP problem. This approach has been used by a number of authors, and a representative example is the control of a continuous stirred tank reactor by Eaton and Rawlings.28 However, NLMPC is computationally demanding, and the optimization problem may not be strictly convex, so that the solution may converge to local minima. If the algorithm does converge to the global minimum, it may take a long time to do so.24 This makes the online implementation of NLMPC a nontrivial task. Instead of NLMPC, the nonlinear quadratic dynamic matrix control (NLQDMC) algorithm29,30 is considered in this work. In NLQDMC, the optimization (described later in this section) is reduced to solving a quadratic program (QP) at every time step. A QP is inherently convex and is thus mathematically more tractable as compared to general nonlinear optimization problems. Consequently, a global solution to the optimization problem can be easily calculated at every time step. The assumption in this approach is that the QP and the original nonlinear optimization problem have similar solutions throughout the batch if the QP is rebuilt at every time step using updated process information. Subsequently, it will be shown in section 5.1 that this is a valid assumption. Because the NLQDMC framework is employed in a shrinking-horizon format, the technique used in this work is termed shrinking-horizon nonlinear QDMC (SNQDMC). 4.1. MPC Design. In this section, the design of an online closed-loop controller based on the MPC philosophy is considered. The essential components of this strategy are summarized below. Objective Function. The 2-norm-squared objective function used in this work is of the following form:

min {||Γy[Y ˆ (k+1|k) - R(k+1|k)]||22 +

∆U(k|k)

||Γu∆U(k|k)||22 + ||Γf[Vf(k) - 20]||22} (13) Vf(k) ) V(k) + P∆U(k|k)

(14)

P ) [p, p - 1, ..., p - (m - 1)]T In eq 13, the first term denotes the projected deviations of the controlled variable, the predicted ethanol trajectory [Y ˆ (k+1|k)], from the time-varying ethanol reference trajectory [R(k+1|k)]; the second term penalizes the manipulated variable moves [∆U (k|k)]; the third term penalizes the deviations of the predicted final volume [Vf(k)] from the upper bound of 20 L. The current volume is given by V(k), whereas the vector P accounts for the sample times for which the input ∆U is applied in order to calculate Vf(k). The relative importance of trajectory tracking, control effort, and meeting the desired final volume can be tuned by altering the values of the diagonal weights, Γy, Γu, and Γf, respectively. Specification of the Reference Trajectory. The optimal input profile that was generated offline in section 3.3 was employed as the ethanol reference trajectory in the closed-loop formulation. Output Trajectory Prediction. The original bioreactor model, discussed in section 2, is a continuous-time nonlinear model. To design a suitable controller for this

process, a controller model is required to predict the process behavior. Analytical linearization of the original nonlinear model was carried out symbolically offline to obtain a function that enabled the calculation of the linearized model online based on the values of the estimated state and past input variables. The local linearization was necessitated by the dynamic nature of the fed-batch operation because the process behavior cannot be represented accurately by linearizing around a single operating point.30 Furthermore, the system is subject to abrupt changes, such as when the original amount of substrate is completely consumed or when the feed is introduced into the reactor, necessitating the reconstruction of locally accurate linear models. The locally linearized models were then discretized at every time step, and step-response coefficients were calculated based on the discretized model. Because the cellularscale reactions are considered to be fast as compared to the reactions at the macroscopic scale, only one characteristic time scale is considered here. This allows the use of a discrete-time step response model to predict the output trajectory. The calculations for discretization and the step response coefficients were performed using the functions c2dmp and ss2step in MATLAB (The MathWorks, Natick, MA). Using these coefficients, the dynamic step response matrix was calculated at every time step to account for the effect of future manipulated variables. The dimension of the dynamic step response coefficient matrix is not a constant, because as the batch proceeds, the prediction horizon p decreases so that the dimension of the step response matrix reduces as well. Considering the input, the controller determines the moves extending until the control horizon, i.e., the next m moves only, so that the vector of current and future manipulated variables is given as

∆U(k|k) ) [∆u(k|k), ∆u(k+1|k), ..., ∆u(k+m-1|k)]T (15) The unmodeled effects at the current sampling time are computed as the difference between the plant measurement obtained from the actual nonlinear model, ym(k), and the model prediction, yˆ (k). In the absence of any further information about these disturbances at future times, it is assumed that the predicted values of the disturbances are equal to the current values over the prediction horizon. With these three contributions, the p-step-ahead predicted output Y ˆ is given as

Y ˆ (k+1|k) ) Y *(k|k-1) + D(k|k) + S uk ∆U(k|k) (16) where Y *(k|k-1) represents the output contribution of past inputs, S uk ∆U(k|k) represents the effect of future manipulated variable moves on the output, and the subscript k denotes that S u is recalculated at every step. The unmodeled disturbance vector of length p is D(k|k). Constraints. Constraints play a very important role in the QP problem solved at every time step of the controller calculation. Because the reactor is constrained to a fixed final volume, the sum of the manipulated variables, i.e., the total substrate feed (plus the initial volume), cannot exceed the final reactor volume bound of 20 L. Additionally, there exists a lower bound constraint that is necessitated by the nonnegativity of the feed into the reactor. Because m future controlled variable moves are considered, the reactor volumes at these future intervals


are calculated as

elements. The manipulated variable rate constraints are then combined into a single convenient expression given as

V(k+1|k) ) V(k|k) + u(k) ∆ts

[ ]

u(k) ) ∆u(k) + u(k-1) so that

Here, the current reactor volume, V(k|k), and the manipulated variable applied at the previous time step, u(k-1), are known at time k and ∆ts is the step size. Similarly, expressions for the future volumes are given as

V(k+2|k) ) V(k|k) + {∆u(k+1) + 2∆u(k) + 2u(k-1)}∆ts

min ∆U(k|k)TH uk ∆U(k|k) - G(k+1|k)T∆U(k|k) (23)

where the Hessian for the QP is T T T u T H uk ) S uT k Γy ΓyS k + Γu Γu + P Γf Γf P

The general expression (when m ) p) for the future volumes is given as

with

[ ] {[ ] [

‚‚‚ 0 ‚‚‚ 0 ... l

] [ ]}

i i - 1 ‚‚‚ 1

∆u(k|k) 1u(k-1) ∆u(k+1|k) 2u(k-1) + l l ∆u(k+m-1|k) iu(k-1)

V(k+i|k) ) V(k|k) + {A∆U(k|k) + B}∆ts

∆ts

∀iep (17)

For m < p, the coefficient matrix is appropriately truncated to have a dimension of p × m. All of these future volumes are constrained by the final reactor volume bound:

V(k+i|k) e 20

∀iep

(18)

This inequality sets a maximum bound for the manipulated variable rate, which in matrix form can be written as

A∆U(k|k) e

20 - V(k|k) -B ∆ts

(19)

Because the flow into the reactor cannot be negative, we have

u(k+j|k) g 0

for j ) 0, 1, ..., m - 1

(20)

u(k+j|k) ) ∆u(k+j) + u(k+j-1) In this manner the lower bound on the manipulated variable rate constraint is expressed in terms of the previous input. In matrix form, this can be written as

IL∆U(k|k) g -U(k-1)

(21)

In this equation, IL is a lower triangular matrix of ones and U(k-1) is a m × 1 matrix with u(k-1) as its

(24)

and the gradient vector is

G(k+1|k) ) 2S

l l

(22)

∆U(k|k)

V(k+3|k) ) V(k|k) + {∆u(k+2) + 2∆u(k+1) + 3∆u(k) + 3u(k-1)}∆ts

V(k+1|k) V(k+2|k) ) V(k|k) + l V(k+i|k)

]

When the prediction equation is substituted into the objective function and the problem is written in the standard QP formulation, the objective function becomes

V(k+1|k) ) V(k|k) + {∆u(k) + u(k-1)}∆ts

1 0 2 1

[

20 - V(k|k) B-A ∆ts IL ∆U(k|k) g -U(k - 1)

with

uT T k Γy ΓyEp(k+1|k)

+ P ΓfTΓf[20 - V(k)] (25)

Ep(k+1|k) ) R(k+1|k) - [Y *(k|k-1) + D(k|k)] (26) Thus, at every sample time, the objective function for the optimization problem is given by eq 23 subject to the constraints given by eq 22. A schematic of the controller design strategy is shown in Figure 5. 5. Simulation Results: Bioreactor Case Study 5.1. Nominal Controller Performance. The SNQDMC algorithm was implemented in MATLAB and the QP was solved online at every controller time step using the function quadprog. The sampling time was 1 h, and the weighting factors used were Γy ) 2, Γu ) 10, and Γf ) 3. The control horizon was set equal to the prediction horizon, i.e., m ) p, and the results are shown in Figure 6. The ethanol concentration at the end of the batch was 38.85 g/L, which is very close to (within 0.8% of) the open-loop final ethanol concentration of 39.17 g/L. The final reactor volume constraint was also satisfied. Upon comparison with the open-loop optimal input profile, it is observed that there is an initial feed addition of 0.13 L in SNQDMC. This is attributed to the final volume penalty term, which contributes to the last terms in the equations for the Hessian and gradient used in the SNQDMC algorithm. Because of the highly nonlinear characteristics of the batch, the final volume predictions using the linear model generated at the initial point are not accurate over the full fed batch. Furthermore, a linearized model is used that is only an approximation of the nonlinear process model. However, as the batch proceeds, better predictions become available and the controller does an admirable job of tracking the reference trajectory. As a comparison, Figure 7 shows the results obtained by solving the nominal problem using NLP. The NLP was implemented in MATLAB, but there are other approaches to solve the NLP problem as well.31,32 Because the focus of the present work is on the control algorithm rather than the efficient solution of the optimization problem, an “off-the-shelf” algorithm was employed. The ethanol concentration at the end of the batch is 38.979 g/L and is greater than that obtained


Figure 5. Schematic of the closed-loop controller design.

Figure 6. Plot of the state trajectories using closed-loop SNQDMC (s) with the closed-loop input profile (-‚-) compared with the openloop optimal trajectory (- -) under nominal conditions (X0 ) 0.2 g/L, S0 ) 100 g/L, V0 ) 1 L, E0 ) 0 g/L, and Si ) 100 g/L). The controller tuning parameters were m ) p, Γy ) 2, Γu ) 10, and Γf ) 3.

from the SNQDMC algorithm (38.85 g/L), which is to be expected. The similarity in the results obtained from the two approaches suggests that the SNQDMC algorithm is a good approximation to the solution obtained from solving the NLP problem. 5.2. Control Horizon Selection. Another interesting feature observed in the simulations was the role played by the control horizon on the final performance, as shown in Figure 8. For a value of m ) 1, the control horizon is limited to the next step. Consequently, the controller calculates only a single move; it is the best single move that the controller can implement in optimizing over the remainder of the batch. This myopic view results in continuous feed addition. As a result, the final ethanol concentration is limited to 33.63 g/L.

As the control horizon length increases, the controller has more input moves available at each time step, and consequently greater flexibility, to optimize over the remainder of the batch, as compared to m ) 1. Thus, the initial feed addition is reduced, and the amount of ethanol obtained at the end of the batch progressively increases. With m ) 3, 36.86 g/L of ethanol is obtained at the end of the batch, whereas for m ) 9, the amount of feed added initially is considerably reduced and results in a final ethanol concentration of 38.84 g/L, which is very close to the 38.85 g/L obtained when m ) p. 5.3. Disturbance Rejection Performance. To test the performance of the controller for disturbance rejection, a +10% change in the glucose feed concentration was made. The sample time and the weighting factors


Figure 7. Plot of the state trajectories using the NLP solution (s) with the corresponding input profile (-‚-) compared with the openloop optimal trajectory (- -) under nominal conditions (X0 ) 0.2 g/L, S0 ) 100 g/L, V0 ) 1 L, E0 ) 0 g/L, and Si ) 100 g/L).

Figure 8. Plot of the feed (left) and volume (right) for different values of the control horizon m ) 1 (s), m ) 3 (‚‚‚), m ) 9 (-‚-), and m ) p (- -) with initial and operating conditions (X0 ) 0.2 g/L, S0 ) 100 g/L, V0 ) 1 L, E0 ) 0 g/L, and Si ) 100 g/L). The controller tuning parameters were Γy ) 2, Γu ) 10, Γf ) 3.

for the objective function were the same as those used in the nominal case, and the results are shown in Figure 9. It is observed that the controller performance deteriorates and attains a final ethanol concentration of only 35.6 g/L. The final volume constraint was also not satisfied. Owing to the increased feed glucose concentration (110 g/L), the volume of feed added in the period from 10 to 16 h, i.e., the initial part of the fed-batch

phase, is reduced. As the batch proceeds, an increased volume of feed is added to account for the volume constraint. Owing to this feed addition at the very end of the batch (17-20 h), a considerable amount of glucose remains unconsumed, which is not an efficient utilization of the available substrate. The final ethanol concentration is reduced as a result of the dilution effect of the added feed volume and the incomplete transfor-


Figure 9. Plot of the state trajectories using closed-loop SNQDMC (s) with the closed-loop input profile (-‚-) compared with the openloop optimal trajectory (- -), for a +10% change in Si. Initial conditions: X0 ) 0.2 g/L, S0 ) 100 g/L, V0 ) 1 L, E0 ) 0 g/L, and Si ) 110 g/L. The controller tuning parameters were m ) p, Γy ) 2, Γu ) 10, and Γf ) 3.

Figure 10. Plot of the state trajectories and feed profile using closed-loop SNQDMC (s) compared with the open-loop optimal trajectory (- -) for a +10% change in Si and re-optimization of the reference trajectory. The controller tuning parameters were m ) p, Γy ) 2, Γu ) 10, and Γf ) 3.

mation of glucose to ethanol. This can be attributed to the fact that the original reference trajectory is no longer optimal and the controller tries to track a trajectory that is not representative of the “disturbed” system. 5.4. Re-optimization Strategy. To address the poor performance in the presence of disturbances, re-optimization can be carried out to generate a new reference trajectory, thereby taking into account the increased

amount of glucose available in the feed. Generally, disturbances occurring in the plant are not known. With the knowledge of the feed added at the start of the batch phase of the reactor run, a system measurement in the batch phase could be utilized to estimate the initial feed concentration that causes the system to deviate from the offline optimal reference trajectory. Because this back-calculation of the feed concentration involves only


Figure 11. Plot of the state trajectories and feed profile using closed-loop SNQDMC (s) with the closed-loop input profile (-‚-) compared with the open-loop optimal trajectory (- -), for a +10% change in Si and a sample time of 0.5 h. The controller tuning parameters were m ) p, Γy ) 4, Γu ) 10, and Γf ) 3.

the numerical integration of the system model, the new (disturbance) feed concentration can be readily obtained. To generate the new reference trajectory, optimization is carried out for the remainder of the batch, from the time point where deviation of the system from the reference trajectory is detected. For the +10% disturbance considered in this work, the system information up to 4 h was utilized. With the knowledge of the feed profile up to 4 h, the glucose concentration in the feed (which in this case is 110 g/L) was calculated as a function evaluation, i.e., the value of Si that satisfies the glucose measurement at 4 h. On the basis of this increased feed glucose concentration, the optimal control policy was recalculated for the remainder of the batch from 4 to 20 h. This new optimal input was used to obtain the modified ethanol reference trajectory, which showed that the new final ethanol concentration target was 43.07 g/L. For the SNQDMC controller, the new ethanol reference trajectory was utilized from a time of 4 h onward. The results of this approach are shown in Figure 10. The controller tracks the new reference trajectory and returns a final ethanol concentration of 42.72 g/L, which is within 0.8% of the re-optimized value of 43.07 g/L. Also, this controller performance was obtained without retuning the controller. The final volume constraint was satisfied, and this closed-loop value is significantly (17%) better than the value of 35.6 g/L obtained with no reoptimization. Furthermore, with the re-optimization strategy, there is a 95% reduction in the deviation of the final ethanol concentration from its optimal value. The re-optimized final ethanol concentration also outperforms the value of 42.19 g/L, which is the final ethanol concentration obtained when the open-loop optimal profile is applied to the disturbed system. 5.5. Effect of Reducing the Sample Time. Ideally, the closed-loop controller should be able to handle the disturbances so that there is no need to re-optimize. The re-optimization strategy also relies on a measurement

that could be noise sensitive. It is observed that the 1 h sample time cannot handle the feed concentration disturbance without resorting to the re-optimization strategy. Thus, it was decided to reduce the sample time by half. This should not be a problem considering the increased computational power that is currently available. Moreover, a 0.5 h sample time would also allow for a much smoother tracking of the reference trajectory than that observed for the 1 h sample time. For the nominal case, the amount of ethanol obtained at the end of the batch was 38.92 g/L with the weighting factors Γy ) 4, Γu ) 10, and Γf ) 3. The results with the disturbance are shown in Figure 11. In this case, no re-optimization was required and the final ethanol concentration obtained was 42.85 g/L. The controller tuning weights were the same as those used for the nominal case. Furthermore, the final volume constraint was also satisfied for both the nominal and disturbance simulations. 6. Summary In this work, the problem of fed-batch fermentation control has been addressed. A model structure that recognizes the multiscale nature of the system was utilized to describe the fermentation of S. cerevisiae on glucose to produce ethanol as a product. The SNQDMC controller design philosophy employed in this work recognizes the dynamic nature of the fed-batch mode of operation, and both the open- and closed-loop control strategies were considered. An open-loop optimal control policy was first determined in order to maximize the ethanol concentration at the end of the batch, with the substrate feed rate as the manipulated variable. In the next step, a nominal controller based on the SNQDMC strategy was designed to track the reference trajectory determined from the open-loop optimization. The controller design also included explicit constraints on both


the state and manipulated variables and demonstrated excellent performance in tracking the reference trajectory while attaining the end-of-batch ethanol concentration. The system performance was shown to be dependent on the prediction and control horizons for the fedbatch reactor operation. The performance of the controller was also evaluated in the presence of a +10% disturbance in the feed concentration, and it was found that online re-optimization to obtain a new reference trajectory significantly improved the disturbance rejection performance of the controller. As an alternative to the re-optimization strategy, a reduced-sample time approach showed improved performance for both the nominal and disturbance cases. Appendix The differences between the model described in the original paper1 and the model used in this work are summarized below:14 (1) The equation for qOmax includes a term with the DOT, so that

(

qOmax ) qˆ Omax 1 +

)(

)

DOT E Ki 1 + DOT

Thus, the maximum respiration capacity is a function of both the ethanol concentration and the amount of dissolved oxygen in the reactor. (2) Initially, qSox is assumed to be equal to qS and qOS is calculated. However, when qOS > qOmax, the flux of qSox is reduced by the factor qOmax/qOS so that qOS ) qOmax. This is an important step to account for the fact that qOS is limited by the maximum respiration capacity qOmax. This leads to the reduction in the fluxes of qSox and qSoxen. (3) The calculation for qEc is performed in an iterative manner, so that

qOmax - qOS E YOE E + KE

qEen )

for i ) 1, i

Closed-Loop Control of Fed-Batch Bioreactors: A Shrinking-Horizon

Recommend Documents