Fiec.espol.edu.ec

Bioprocess Modelling for Learning Model
Predictive Control (L-MPC)
Mar´ıa Antonieta Alvarez1,, Stuart M. Stocks2, and S. Bay Jørgensen1, 1 CAPEC, Department of Chemical and Biochemical Engineering, Technical University of Denmark, Søltofts Plads, 2800 Kgs. Lyngby, Denmark 2 Novozymes, Bagsværd, Denmark Abstract. Batch and Fed-Batch cultivation processes are used exten-
sively in many industries where a major issue today is to reduce the
production losses due to sensitivity to disturbances occurring between
batches and within batches. In order to ensure consistent product quality
by eliminating the influence of process disturbances it is very important
to consider implementation of monitoring and control and thereby signif-
icantly improve the economic impact for these industries. A data driven
modeling methodology is described for batch and fed batch processes
which is based upon data obtained from operating processes. The chap-
ter illustrates how additional production experiments may be designed to
improve model quality for control. The chapter also describes how the de-
veloped models may be used for process monitoring, for ensuring process
reproducibility through control and for optimizing process performance
by enforcing learning from previous batch runs through Learning Model
Predictive Control (L-MPC).
Fermentation processes are spreading into traditional chemical industries as moreenvironmentally friendly production methods gain importance. Most often fer-mentations are carried out as batch or fed-batch cultivations which have beenused extensively in bio-technical and pharmaceutical industries. However micro-bial cultivations exhibit significant sensitivity to disturbances occurring betweenbatches and within batches which as evidenced by a relatively large batch tobatch variability often experienced in production. Therefore there is significantinterest in development of methods which may ensure more consistent productquality by eliminating or counteracting the influence of process disturbances. Inorder to ensure more consistent product quality it is very important to considerimplementation of monitoring and control, and thereby significantly improve theeconomic impact for these industries.
 Present address: Facultad de Ingeniera en Electricidad y Computaci´on. Escuela Su- ecnica del Litoral. Guayaquil. Ecuador.
 Corresponding author.
M. do Carmo Nicoletti, L.C. Jain (Eds.): Comp. Intel. Tech. for Bio. Mod., SCI 218, pp. springerlink.com  Springer-Verlag Berlin Heidelberg 2009 M.A. Alvarez, S.M. Stocks, and S.B. Jørgensen Given the main goal of reducing batch to batch variability the question then is which possibilities there are or may be developed for satisfying this purpose.
When viewed from the properties of microorganisms then these cultivations areinfluenced by many variables which render the control problem multivariatefrom the outset. Furthermore a control design which enables tracking the de-sired batch trajectory within some margin is highly desirable in order to avoidthat the cultivation triggers a one of the many stress responses possible andthereby deviate significantly from the desired behavior. Hence a multivariatecontrol design which enables control actions relatively often throughout the fed-batch production phase would be desirable. However multivariate control callsinevitably for a model based control design. Thus the question is how to developa model which can be suitable for such a control design?. Given the presentknowledge of microorganisms and especially the limited understanding of thebehavior of their regulatory network precludes a first principles modeling ap-proach. Hence a data driven approach seems to be the only viable alternative.
Therefore such an approach is investigated in this chapter.
The batch processing types covered in this chapter includes Batch, Fed-batch and periodic operation which all have the common traits of a repeated opera-tion which start from nearly the same initial conditions. Thus the time withinthe batch and the batch number are the two characteristic independent vari-ables. Batch processing is subject to variations in raw material properties, instart-up initialization and other disturbances during execution. These differentdisturbances introduce often significant variations in the final product quality.
Compensating for these disturbances have been difficult in the past due to thenonlinear and time-varying behavior of batch processing and to the fact thatreliable on- or in-line sensors for monitoring final product quality rarely areavailable. Consequently development of a systematic methodology which canensure reliable reproducible operation may provide significant benefits for batchprocessing.
Each batch operation may be defined as a series of operational tasks, i.e.
mixing, reaction and separation. Within each task a set of subtasks, e.g. heat-ing/cooling, (dis-)charging is handled. There may be more than one feasible setof operational tasks that can produce the specified product(-s). Consequentlyan optimal sequence of tasks and subtasks with respect to a defined objectiveneeds to be identified. This set of operational tasks is labeled the optimal batchoperations model. Thus the Batch Operations Model combines the batch pro-cessing tasks normally specified in a generic recipe with the batch equipmentunder availability and other resource constraints.
The question of which methodology to select is based modeling using general empirical models combined with the possibility for executing control throughoutthe fed-batch phase. A method based on prediction of end of batch propertiesconsiders midterm adjustment, which may be too late for many industrialcultivations. Therefore experimental adaptation (optimization) of the Batch Op-erations Model, e.g. the method based on "the solution model" through trackingconstraints such as the Necessary Conditions of Optimality or Bioprocess Modelling for Learning Model Predictive Control (L-MPC) determination of the presently (most) constrained variable may be feasible.
However for the present study the Grid of Linear Model method wasselected since that method actually is quite general in its approach.
This chapter presents methodologies based upon the previous ideas for batch or periodic operation to provide a predictive modeling which can enable reliablereproducible operation using control and which may enable optimizing oper-ation. The contribution comprises data driven time series modeling of batchprocesses based upon a large set of locally linear models and a learning modelpredictive control methodology. The modeling methodology produces both aLinear Time-Invariant state space model representation for inter-batchprediction and a Linear Time-Varying state space model representationfor intra-batch prediction. The modeling approximates the non-stationary andnonlinear behavior of batch processes with a set of local but interdependent lin-ear regression models parameterized as AutoRegressive Moving Average modelswith eXogenous inputs. Tikhonov Regularization is applied toreliably estimate the many parameters of this model set. However since the modelis intended for control it is important to ensure that the model is suitable forthis purpose. Therefore design of identification experiments with input pertur-bations is proposed in this chapter and investigated both in simulation and on anindustrial cultivation. The resulting model may be applied for Learning ModelPredictive Control . This control methodology is presented for controlof repeated operation of stochastic Linear Time-Varying systems with finite timehorizons together with tuning requirements for ensuring guaranteed convergenceand hence closed-loop stability. The methodologies have been implemented as aMatlab toolbox named Grid of Linear Model . The application of thesemethods is illustrated through both the simulation study and the experimentalpilot plant investigation.
The chapter is structured such that the methods are briefly introduced and subsequently they are illustrated on a simulation benchmark example. Here theexperimental design, the model estimation and the design and perfor-mance. Subsequently an industrial cultivation is modeled with the aim to imple-ment an control design. Here both the model estimation and implemen-tation of a Kalman filter for estimating unmeasured variables, especially enzymeactivity are illustrated. The implementation of a control design is brieflydiscussed, before the conclusions are drawn.
This section briefly presents the different methods used in this chapter. For mod-eling estimation of a fermentation process, the methodology applied is Grid ofLinear Model as explained in subsection The Learning Model Pre-dictive Control methodology used for the batch fermentation processis presented in subsection To develop models from a data set consist-ing of a relatively low number of batches it is desirable to design experimentswith perturbations in the actuator variables. The experimental design methodapplied is presented in subsection M.A. Alvarez, S.M. Stocks, and S.B. Jørgensen Modeling for Control
Batch processes are modeled with the toolbox as a sets of N models.
Such a set of models could also be referred to as one batch model. Thesemodels can be parameterized in a number of ways, but in the present contri-bution an parameterizations was chosen. This choice of parametriza-tion offers a simple multivariable system description with a moderate numberof model parameters. The objective of the model set is to quantify the causalcorrelations between the process outputs yk,i ∈ IRny , inputs uk,i−1 ∈ IRnu, anddisturbances vi ∈ IRny , for i = 1, . . , t, at times t = 1, . . , N in batch k. To sim-plify notation, define the input uk, output yk, shifted output y0, and disturbance vk profiles in batch k as uk = uk,0 uk,1 . . uk,N−1 yk = yk,1 yk,2 . . yk,N y0 = y k,0 yk,1 . . yk,N −1 vk = vk,1 vk,2 . . vk,N The toolbox models the differences between two successive batches Δyk = yk − yk−1 = −AΔy0 + BΔuk − Cvk where Δuk = uk − uk−1 and Δy0 = . The batch model may be converted into different representations dependent on the particularapplication task. If the task at hand is to predict (or simulate) the behavior ofa batch before it is started, the Eq. is convenient Δyk = HΔyk,0 − GΔuk + F vk where Δyk,0 = yk,0 − yk−1,0. The change in the initial conditions Δyk,0 can beconsidered as either an input/control variable or a disturbance. Eq. is alsoconvenient for the task of classification (e.g. normal or not) of a batch after ithas been completed. Furthermore, Eq. can be used to determine open-loopoptimal recipes in the sense of optimizing an objective for the batch. If the objective is to minimize the deviations ek = ek,1, ek,2, . . , ek,N , ek,t ∈ IRny,from a desired Batch Operations Model ¯ y, then Eq. can be modified into y − yk = ek−1 − HΔyk,0 + GΔuk − F vk The two Eqs. and of the batch model above are applicable to off-line or inter-batch type applications. For on-line estimation, monitoring,feedback control, and optimization however, it is convenient to use a state spacerealization of the batch model.
In an observer canonical form, which is structurally a minimal realization, the state space realization is given as xk,t = Atxk,t−1 + BtΔuk,t−1 + Etvk,t Δyk,t = yk,t − yk−1,t = Cxk,t Bioprocess Modelling for Learning Model Predictive Control (L-MPC) with Δuk,t−1 = uk,t−1 − uk−1,t−1 and the state space model dimension nx =ny max (ni(t) 1 ≤ t ≤ N, i = A, B, C) and the initial condition xk,0 = CΔyk,0.
The state space model matrices At ∈ IRnx×nx , Bt ∈ IRnx×nu(t), Et ∈ IRnx×ny ,and C ∈ IRny×nx contain the corresponding block columns in the batch ARMAXmodel (A, B, C). To exemplify, assume that nA(t) = nB(t) = nC(t) + 1 = 3 andthe state space model matrices will be given as −at,t−1 I 0 −at+1,t−1 0 I , Et = −ct+1,t −at+2,t−1 0 0 C = I 0 0 Just as Eq. the state space model form is convenient for prediction, monitoring, and optimization type applications, and it facilitates on-line imple-mentation of such applications. Furthermore, the state space model form is particularly well suited for closed-loop or feedback control applications. Fortracking control applications the state space model form can be modified to xk,t = Atxk,t−1 + BtΔuk,t−1 + Etvk,t yt − yk,t = ek−1,t − Cxk,t In order to use Eq. for tracking control, it is necessary to estimate the states x) based on noisy observations of the outputs. Assume that during a batch (k), observations zk,t of the outputs yk,t are collected at times t = 0, 1, . . , N andlet the optimal estimate of the state xk,t in batch k at time t1 given data up to and including time t2 be given as the conditional mean (e.g. obtainable with aKalman Filter) ˆ } where the information Ik,t = {zk,t, Δuk,t, Ik,t−1}, Ik,−1 = Ik−1,N , I0,−1 = {y−1, u−1, z−1} (7) Then the tracking error ek,t in batch k at time t1 given data up to and includ- ing time t2 is estimated as ˆ where the smoothed ek−1,t1 N − C ˆxk,t1 t2 estimate of the error profile in batch k − 1 is given as k−1 k−1 = k−1,N N yk−1 k−1 is the smoothed output profile estimate from batch k − 1 (e.g.
obtainable with a Kernel Smoother). Note that the superscript signifies anestimate.
Learning – Model Predictive Control
Given the previous developed models for batch operation then it is possibleto develop a control paradigm which handles both intra batch and inter batchdisturbances. This is achieved by combining Model Predictive Control, which M.A. Alvarez, S.M. Stocks, and S.B. Jørgensen handle intra batch disturbances, with a Learning Control, which handle batchto batch disturbances. To utilize available information during a batch to obtainthe best possible tracking performance, the following Learning Model PredictiveControl formulation k,l,t}N −1 = arg min ek,i t + Δu k,i }N −1 i=t+1 k,i t xk,i t = Ai ˆ xk,i−1 t + BiΔuk,i−1 ek,i t = ˆek−1,i N − C ˆ k,i−1 + uk−1,i−1 ≤ umaxi−1 ymin ≤ ¯y ek,i t ≤ ymaxi is solved at times t = 0, 1, . . , N − 1 in batch k and the control sequence Δuk = Δuk,0,0 Δuk,1,1 . . Δuk,N−1,N−1 then approximates a closed-loop optimal control sequence. I.e., at time t inbatch k, Eq. is solved based on updated state estimates, and the inputuk,t = Δuk,t + uk−1,t is implemented on the process.
Experimental Design for Control
In order to efficiently reject disturbances in bioprocesses it is important to be ableto act as soon as deviations from a desired performance is detected, otherwise thebioprocess may develop into operating regions from which it is not possible torecover. One such situation may occur if the microorganisms experience a lack ofoxygen, which may lead to excessive cell death. Therefore it is essential to be ableto apply control throughout the batch operation to enable tracking of the desiredoperations model. However this desire implies that the data driven modelingmust be based upon data where the actuator variables have been moved duringthe process to provide sufficient information for model parameter estimation.
Thus experimental design for identification becomes very important for providingsuch informative data.
The subject of experimental design is well known and studied for continuous processes, e.g. , however, for batch processes additional aspects come intoplay, just as for model estimation as treated briefly above. These aspects wereinvestigated by The key point is that by introducing perturbations aroundthe nominal Batch Operations Model it is possible to obtain informative datafor identification of batch models. These include a.o.: 1. The batch operation should be perturbed around the nominal batch opera- tions model for the potential control actuator variables.
2. The magnitude of the perturbations should ensure a balance between prox- imity to the nominal batch and sufficient perturbation to obtain informativedata, as judged, e.g. by the Fisher information matrix.
3. For noise-corrupted measurement data, the noise content can be reduced by smoothing, e.g. Kernel Smoothing.
Bioprocess Modelling for Learning Model Predictive Control (L-MPC) The detailed design of the perturbations is treated in connection with the benchmark case and the industrial pilot plant case, described in Section andSection respectively.
The methodology is used to develop a discrete-time state space model rep-resenting a batch process based on historical operating data. In subsection a simple Benchmark case is used to represent a fermentation process (PenicillinFermentation). From the simulation of the benchmark model a set of batches isobtained, with this data set a state space model is estimated using the methodology (subsection In subsection the methodology is ex-plained for control design. Based upon the obtained model a Learning ModelPredictive Control controller is designed (subsection and demon-strated through simulation under two scenarios: one with inter batch distur-bances and another with intra batch disturbances. This simple benchmark caseis used to illustrate the methodology behind and to investigate modeldevelopment for prediction of process behavior using a data set with a limitednumber of batches for modeling. In subsection the conclusions of the bench-mark case are given.
Benchmark Simulation Model
The example fed batch process is simulated penicillin fermentation. The modelis a modified version of the model in which describes the growth of biomassand formation of the single product from a pure substrate simulated with theparameters from , these values are shown in Table The model equationsare: dX = μX − FX F (S = − μX − θX − M dP = θX − KP − FP dV = F for t ∈ [t0, tf ], where X (g/L) is the biomass concentration, S (g/L) is the sub-strate concentration, P (g/L) is the product concentration, V (L) is the reactorvolume, F (L/h) is the feed flow rate, YX and YP are yield coefficients andSf (g/L) is the substrate feed concentration. MX represents a constant specificmaintenance demand of the cells and K represents a constant first-order de-cay rate for the product. For the kinetic equations, specific growth rate μ uses M.A. Alvarez, S.M. Stocks, and S.B. Jørgensen Table 1. Simulation model parameters
K21 0.0001 gl MX 0.029 h−1 X0 0.01 h−1 S0 μmax 0.11 h−1 P0 θmax 0.004 h−1 V0 Michaelis Menten kinetics while the product formation rate θ were modeled withContois kinetics : μ = μmax S + K1X θ = θmax K22S2 + S + K21 This simple model (physical model Eqs. - provides the example process used to generate a data set for model identification using the methodology(black box model). After the black box model is development, the same physicalmodel is used to simulated the "real process" under control.
The penicillin simulation case represents a cultivation with two input variables (substrate feed concentration Sf and feed flow rate F ) and four measured outputvariables (biomass X, substrate S, product P (penicillin) and reactor volume V ).
The simulation generates a data set with 7 batches where every batch runs for150 hours with measurements sampled once per hour (150 samples). The numberof batches was selected as a realistic case to test the estimation capability of themethodology to a realistic case, where perturbations are used to ensurethe information content in data. The Penicillin simulation starts with a batchphase in which substrate, nutrients and inoculum have been loaded into thesterilized reactor where only oxygen is supplied during cultivation. Initially in thebatch phase the nominal value of substrate concentration is high (80 g/L) whilethe biomass formation rate reaches its highest value after an almost exponentialgrowth. The initial substrate composition and amount varies with a standarddeviation of 20% between the seven batches. When the substrate concentrationfalls below 5 g/L the fed-batch phase is set to start, the feed flow rate is 2 L/h andthe substrate feed concentration is 400 g/L, the initial reactor volume is 250 L.
Independent Pseudo Random Binary Signal sequences are used with a clock period of 3 hours to perturb each input variable. The perturbation am-plitude for the feed flow rate is 1 L/h, whereas, the perturbation amplitude forthe substrate concentration in the feed is 200 g/L. The amplitude of both sig-nals is chosen sufficiently large to generate measurable process variations whichare sufficient to obtain reliable information on the process behavior but not solarge that the process behavior drifts away from its normal behavior. Figure shows the input sequences to the process with perturbation for the seven batchessimulated.
Bioprocess Modelling for Learning Model Predictive Control (L-MPC) Substrate feed concentration Fig. 1. Feed flow rate and Substrate concentration signals with disturbance. Top: Feed
flow rate signal vs time; bottom: Substrate concentration in the feed flow vs time.
The outputs of the process for theses inputs are depicted in Figure It is clearly possible to see the effect which the perturbations generate in the processvariables.
To investigate the influence of the perturbation amplitude for model estima- tion, three data sets again each with seven batches have been simulated with dif-ferent perturbation amplitudes. The perturbations that characterize each datasetare described in Table Table 2. Characterization of the perturbations for each dataset in the seven simulated
benchmark experiments
Amplitude perturbed % of the Signal Feed flow rate (L/h) Substrate feed concentration (g/L) Benchmark Model Estimation and Validation
The following sections are based on . A model is estimated using a novel in-terpretation of Tikhonov Regularization, the estimated model is validated onindependently perturbed simulations. Due to the repetitive nature and the finitehorizon, the measurements are collected from a grid of sample points in time.
With these samples points it is possible to model the evolution between two con-secutive sample points in a batch with local models, which are labeled grid-point models (Figure Estimation of the parameters in the grid point models M.A. Alvarez, S.M. Stocks, and S.B. Jørgensen Fig. 2. Results of the simulation. Top left: Biomass conc. vs time; top right: substrate
conc. vs time; bottom lef: product conc. vs time; bottom right: volume vs time.
is possible because at each sample time the measurements are repeated as manytimes as there are batches. Clearly however the Linear Time-Invariant model which represent the behavior from batch to batch is rather large since thenumber of parameters is proportional to the product of the number of measure-ment intervals multiplied with the number measured process variables plus thenumber of input variables. Therefore it is essential to use a sound regularizationmethod which ensures that only parameters which extract information from themeasured data are actually estimated. This is achieved by using regularization.
Note that when considering the behavior within a batch then the batch processmay be represented as a Linear Time-Varying batch model.
Input-Output Modeling: The model parametrization may be carried out in
two ways, using either AutoRegressive models with eXogenous inputs or
AutoRegressive Moving Average models with eXogenous inputs. The
software gives the option to use different parameter estimation methods,
i.e. as Least Squares Ridge Regression and Tikhonov Regularization
. Here is preferred, since it provides a stable and efficient method for
determining the model parameter estimates.
To illustrate the modeling concept the different formulations of the and models used for modeling is presented next.
Bioprocess Modelling for Learning Model Predictive Control (L-MPC) Fig. 3. The three dimensions of batch data for one particular variable, e.g. biomass
concentration. The filled circles represent the sample points taken every 10 hours. The
smooth connecting lines symbolize the development of the particular variable through
the batch. The connecting lines between batches illustrate the development form batch
to batch. Note the finite duration of the batch.
The inputs and outputs of the batch operation are defined as deviations from a reference trajectory for the batch. The reference trajectory is indicated withan overbar. Thus the model at each discrete time instant t is given as: – Input variable ut ∈ IRnu(t) with reference ¯
ut ∈ IRnu(t) – Output variable yt ∈ IRny(t) with reference ¯
yt ∈ IRny(t) – Disturbance variable wt ∈ IRny(t)
Using an model parametrization, the output deviation from the reference yt − yt at time t may be given in the following locally linear equation yt − yt = −at,t−1(¯ yt−1 − yt−1) − . . − at,t−nA(t)(¯yt−nA(t) − yt−nA(t)) ut−1 − ut−1) + . . + bt,t−n B (t)(¯ where nA(t), nB(t) ∈ [1,. .,t] are the structured grid-point model ordersand ai,j ∈ IRny(i)xny(j) and bi,j ∈ IRRny(i)xnu(j) are the structure grid pointmodel parameter matrices, wt is the disturbance. The model type or , depends on how the disturbance is modeled. For the entire batchthe input-output model is expressed as y − y = −A(¯ y0 − y0) + B(¯ u − u) + w M.A. Alvarez, S.M. Stocks, and S.B. Jørgensen where the input u, output y, shifted output y0, and the disturbance profile ware u = [u y = [y y0 = [y w = [w The process disturbances are caused by several influences which include bias in the reference input profile ¯ u, the effect of process upsets, modeling errors from linear approximations and errors due to bias in transition times between sets ofactive constraints. The disturbance profile w contains contributions from bothbatch wise persistent disturbances and random disturbances with no batch wisecorrelation. Hence, the disturbance is modeled as a random walk model withrespect to the batch index k wk = wk−1 + Δwk The difference between two successive batches is given by Δyk = yk − yk−1 = −A(y0 − B(uk − uk−1) − Δwk the disturbance profile increment Δwk is modeled with a Moving Average model with respect to time Δwk,t = vk,t + ct,t−1vk,t−t + · · · + ct,t−nC(t)vk,t−nC(t) where the model order nC(t) ∈ [0, . . , t − 1]. In matrix form the disturbancemodel is: Δwk = Cvk Then the following model can describe the difference profile between two successive batches Δyk = −A(y0 − y0k−1) + B(uk − uk−1) − Cvk A key difference between and models is that an model includes the model of the disturbance. Depending on the particular applicationpurpose the model (Eq. can be converted into different representa-tions. The following Eq. is convenient if the task is to simulate the behaviorof a batch or for classification (e.g. normal or not) of a batch after it has beencompleted, an other task can be to determine open-loop optimal recipes in thesense of optimizing an objective for the batch.
Δyk = HΔyk,0 − GΔuk + F vk The form of the batch model (Eq. is used for off-line or inter- batch type applications, while for on-line application such as on-line estimation,monitoring, feedback control and optimization; it is convenient to use a SpaceState realization of the batch model: xk,t = Atxk,t−1 + BtΔuk,t−1 + Etvk,t Δyk,t = Cxk,t Bioprocess Modelling for Learning Model Predictive Control (L-MPC) with the initial condition xk,0 = CΔyk,0 This Space State model is precisely the form given in Eq. after which the ma- trices are defined. Like Eq. the Space State model form (Eq. is convenientfor prediction, monitoring and optimization type applications. In addition theSpace State model form (Eq. is well suited for close-loop control applications.
Parameter Estimation: Tikhonov Regularization provides an efficient
method to obtain a trade off between bias and variance of the model parameter
estimates. It is a coefficient shrinkage based parameter estimation method that
can incorporate model properties into weighted Least Squares estimates.
The formulation for the estimation problem used in is
V , W, Λ) = arg minθ JTR
s.t. JT R = 
+ L
Y − Xθ2W where L is a structured penalty matrix which maps the parameter vector θ into
the desired parameter differences, Λ is a diagonal weighting matrix that weights
the parameter differences. The penalty matrix L consists of five sub-matrices
L = [L
L ]
the five sub-matrices are weighted individually by block diagonal weightingmatrices Each of the five sub-matrices Li for i = 1, ., 5 imposes a penalty on certain
properties of the parameters as briefly described in the following.
L1: Penalizes the model parameter time evolution by penalizing the approximate
first order time derivative of the parameters θ. This means that the quantities
(using A matrix as an example)
at,t−l(i, j) − ˆat+1,t−l+1(i, j) for i, j = 1, ., ny(t), t = 1, ., N − 1 and l = 1, ., nA(i, j, t) are penalized ac-cording to the corresponding scalar weight λa L2: Penalizes non-smoothness of the model parameter time evolution by penal-
izing the approximate second order time derivative of the parameters θ. This
means that the quantities (using A matrix as an example)
at−1,t−l−1(i, j) − 2ˆat,t−l(i, j) + ˆat+1,t−l+1(i, j) M.A. Alvarez, S.M. Stocks, and S.B. Jørgensen for i, j = 1, ., ny(t), t = 2, ., N − 1 and l = 1, ., nA(i, j, t) − 1 are penalizedaccording to the corresponding scalar weight λa L3: Penalizes the time evolution of the impulse responses by penalizing the ap-
proximate first order time derivative of the impulse response of the local models
θt. This means that the quantities (using B matrix as an example)
bt+l−1,t−1(i, j) − ˆbt+l,t−1(i, j) for i = 1, ., ny(t), j = 1, ., nu(t), t = 1, ., N − 1 and l ∈ { l nB(i, j, t + l) >l, l ∈ [1, ., N − t]} are penalized according to the corresponding scalar weightλb L4: Penalizes the non-smoothness of the impulse responses by penalizing the
approximate second order time derivative of the impulse response of the local
models θt. This means that the quantities (using B matrix as an example)
bt+l−1,t−1(i, j) − 2ˆbt+l,t−1(i, j) + ˆbt+l+1,t−1(i, j) for i = 1, ., ny(t), j = 1, ., nu(t), t = 1, ., N − 2 and l ∈ { l nB(i, j, t + l) >l, l ∈ [1, ., N − t]} are penalized according to the corresponding scalar weightλb L5: Penalizes the variance of the model parameter estimates ˆ
θ (L5 = I). This
means that the quantities (using C matrix as an example)
for i, j = 1, ., ny(t), t = 2, ., N and l = 1, ., nC(i, j, t) are penalized accordingto the corresponding scalar weight λc . It is noted that if L = L
problem (Eq. is reduced to a standard form problem or a problem.
As for other regularization methods such as Ridge Regression, this method introduces bias to the model parameter estimates. The trade off between variance
and bias of the parameter estimates determines the predictive capabilities of the
model. The weighting matrix Λ determines the coefficient shrinkage. Selection of
which elements to include and to what extend constitutes important parameters
for the design of an optimal model for a given purpose.
Model Selection: The quality of the model is evaluated by how well the model
generalizes to an independent data set. The method used is cross-validation,
which estimates the generalization error G when a model ˆ
θ is applied to an independent data set .
k − Δˆ where Δˆ yk is the model prediction of the regressors in batch k of the independent validation data set and ¯ W is a block diagonal weighting matrix with symmetrical, Bioprocess Modelling for Learning Model Predictive Control (L-MPC) positive definite, block elements ¯ Wt ∈ IRny(t),ny(t) and the weighting matrix ¯ should be chosen as Wt(i, i) = ((N − 1)nymax { Δˆ yk,j(i) k = 1, ., NB; j = 1, ., N })−1 (37) for i=1,.,ny and t=1,.,N. The size of the validation data set N val is typically 50% of the size of the data set used for model parameter estimation N est.
In practice the validation is performed through a combination of visual in- spection of pure simulation, i.e. prediction of the entire batch. This is supportedwith comparison of model performance using a scalar measure of the FIT, orrather the misfit, which also is used by for quantitative analysis of thequality of the model on validation data. Thereafter the model to be used for thespecific modelling purpose is selected.
Estimated Benchmark Model: As mentioned, there are three data sets (a, b
and c - see Table 1)each containing seven batches to be modeled; the objective of
the modeling is to determine the influence of the input perturbation amplitude
on model identification when only a few batches are available. In the user
selects the model type or the maximum model order and the
estimation method (weighted Least Squares , Ridge Regression or
Tikhonov Regularization ).
To model each of the three data sets, splits up the data into an esti- mation set, a model validation set and a test set. Batches 1-3 have been used foridentification, batches 4-5 have been used for validation and batches 6-7 havebeen used for model test. The model type used is ARX with a maximum modelorder of two and the estimation method used is can estimate themodel using Pure Simulation and One Step Ahead prediction er-ror profiles: estimates the outputs based on initial conditions at the start ofthe cultivation, this means that an early error will remain throughout the batch.
Whereas in estimation, the output estimated is based on using the measuredoutput as initial condition. To provide a compromise between these two types ofprediction, then another user parameter is introduced Model Purpose whichis a number between 0 and 1, where 0 estimates the model using One StepAhead Prediction while 1 is for Pure Simulation Prediction. Thus the ModelPurpose weights between prediction of a one-step ahead error and a puresimulation - where the whole batch is predicted - error. The larger the valuethe closer the estimate is to a pure simulation prediction. The default setting inSoftware is the golden cut = 0.618.
The models estimated will be compared using a Pure Simulation Prediction generalization error (Eq. as calculated by the Software. This is com-plemented with visual inspection of the model parameters and predictions. ThisFIT (or rather misfit) measure is between 0 and 1, where near 0 means that themodel predicts the process trajectory very well, whereas a FIT near 1 impliesa poor prediction of the process behavior. The parameters for the estimationsetting in are shown in Table for the three data sets (a, b and c - seeTable .
M.A. Alvarez, S.M. Stocks, and S.B. Jørgensen Table 3. Parameters for GoLM estimation for the three data set (a, b and c).
Model Structure Max. model order Model Purpose Tikhonov Regularization Table 4. FIT for Pure Simulation Prediction for ARX models
Data set FIT Biomass FIT Substrate FIT Product FIT Volume FIT Total With the FIT of the three models (see Table it is possible to conclude that when there are few batches available to be used for modeling, it is extremely im-portant that there is sufficient information in the data for identification. Thus,it is essential to realize designed experiments for model identification. With per-turbation in the inputs, the outputs contain information about the dynamicsfrom these inputs and such data is sufficiently good to be used for modeling forcontrol purposes.
The criteria used to choose the model is the FIT obtained by each model as listed in Table and visual inspection of Figure For the first data set (modela) (Figure a)), the perturbed amplitude is 10% of nominal value of the signal.
The model obtained does not simulate the real process well, hence the perturbedinputs do not generate sufficient information in the outputs. Whereas for dataset model b (Figure the perturbed amplitude is 25% of the nominal valueand with this perturbation in the inputs, it is possible to obtain a good modelof the process. The data set model c -with a perturbed amplitude of 37% of thenominal value- in general, only provides a little improvement in the total FIT.
But Figure shows that the model estimated for substrate develops negativevalues during the batch phase. Thus, perturbation amplitude for identificationis 25% of the nominal value (data set b) provides the most satisfactory data formodel identification of this benchmark case.
The graphics that presents Software for the estimation results show the comparison between measurement, Pure Simulation Prediction, One Step AheadPrediction and the previous batch. Note that a naive prediction represents theprevious batch. Prediction represents the most demanding application of abatch process model. Hence, if a model quality based on model predictionis reasonable in terms of the estimated generalization error, then the modelwill have credibility in any similar application. There is a tendency that modelselection based on Prediction under-smoothes, while Prediction errorhas a tendency to over-smoothing. As Figure show, for model b providesclearly the most reasonable prediction. Obviously OSA provides a seeminglybetter FIT, however for such a model the prediction horizon is basically limitedto just one sample.
Bioprocess Modelling for Learning Model Predictive Control (L-MPC) Pure Simulation Prediction (a) model a
Pure Simulation Prediction (b) model b
Pure Simulation Prediction (c) model c
Fig. 4. (a): Simulation of batch 6 using a model obtained from the model estimation
using the data set a; (b): Simulation of batch 6 using a model obtained from the model
estimation using the data set b; (c): Simulation of batch 6 using a model obtained from
the model estimation using the data set c
M.A. Alvarez, S.M. Stocks, and S.B. Jørgensen For the previous model estimation (model estimated a, b and c). The analysis was done keeping the = 0.618 (see Table and varying the input perturbationamplitude for the three data sets (a, b and c - see Table The model obtainedwith data set b (model b) is attempted further improved using the decreasingweight of Pure Simulation as shown in Table Table 5. Parameters for model estimation for data set b for estimating three
models with different values of Model Purpose
Model name Model Structure Max. Model order Model Purpose Tikhonov Regularization Tikhonov Regularization Tikhonov Regularization Figure visualizes the obtained FIT for the three models (model 1, 2 and 3), while Table presents the FITs for each variable in each model. Clearlymodel 1 in Figure (a) shows a significant improvement of the model FIT. Thisinvestigation illustrates that the Model Purpose parameter can be used toimprove the estimated model according to the model application which here isto predict the entire batch duration.
Table 6. FIT for Pure Simulation Prediction of batches 6-7 for ARX model using data
set b with different value of
Model name FIT Biomass FIT Substrate FIT Product FIT Volume FIT Total Batch Control Design
This subsection introduces the use of Learning Model Predictive Controlin batch processes, first (subsection introduces of ap-plying modeling. Then, subsection presents an implementation of thecontroller on the Benchmark example.
The idea behind Learning Model Predictive Control is tocombine the asymptotic performance convergence of Iterative Learning Controlwith the approximately closed-loop performance of Model PredictiveControl . proposes a Learning Model Predictive Control of batch pro-cess which is described by stochastic Linear Time-Varying system directlyobtained from Eq. : xk,t = Atxk,t−1 + BtΔuk,t−1 + Etvk,t yt − yk,t = ek−1 − Cxk,t Bioprocess Modelling for Learning Model Predictive Control (L-MPC) Pure Simulation Prediction (a) model 1
Pure Simulation Prediction (b) model 2
Pure Simulation Prediction (c) model 3
Fig. 5. (a): Simulation of batch 6 using model 1; (b): Simulation of batch 6 using model
2; (c): Simulation of batch 6 using model 3
M.A. Alvarez, S.M. Stocks, and S.B. Jørgensen for t = 1, . . , N with initial condition xk,0 = −Cvk,0 This means that the state xk,t ∈ IRnx of the system at time t in batch k, is given by linear mappings of the state xk,t−1 at time t − 1 in batch k; thecontrol correction Δuk,t−1 = uk,t−1 − uk−1,t−1 ∈ IRnu(t−1) and a zero-meanGaussian disturbance vk,t ∈ IRny . The tracking error ek,t ∈ IRny of the systemis given as the difference between the system output yk,t ∈ IRny and thedesired output reference y¯k,t ∈ IRny .
The formulation for the is k,l,−1}N −1 = arg min ek,i t + Δu k,i }N −1 i=1 k,i t xk,i −1 = Ai ˆ xk,i−1 −1 + BiΔuk,i−1 ek,i −1 = ˆek−1,i N − C ˆ Δuk,i−1 + uk−1,i−1 ≤ umax i−1 ymin ≤ ¯ yi − ˆek,i −1 ≤ ymax i if the weighting matrices Qk and Rk are block diagonal Qk = diag(Qk,t) Rk = diag(Rk,t) and the weighting matrices Qk,t and Rk,t are all symmetric and positive definite.
proves that the optimal solution Δuk,−1 = {Δuk,l,−1}N−1 to the satisfies the stochastic design requirement and guarantees convergence (seefor details).
Thus, the cost function J defined is (Eq. Where t = −1, 0, . . , N and −1 ≤ j < N.
Jk(ˆek,t(Δuk,j), Δuk,j) = ˆek,t(Δuk,j)2 + Δu State and Output Estimation: For state and output estimations, implements the Kalman Filter recursive equations
xk,t t = ˆ xk,t t−1 + Kk,t(zk,t − C ˆ xk,t t−1 − ˆ yk−1,t N ) xk,t t−1 = AS ˆ t xk,t−1 t−1 + BtΔuk,t−1 + EtStR−1 k,t−1 − ˆ yk−1,t−1 N ) for t = 0, 1, . . , N , with initial condition xk,0 −1 = 0 The Kalman Filter gain matrix Kk,t and state estimate covariance matrix Pk,t tare propagated by the recursions = Pk,t t−1(I − CK ) Kk,t = Pk,t t−1C(CPk,t t−1C + R + ¯ k−1,t N )−1 Pk,t t−1 = AStPk,t−1 t−1AS Bioprocess Modelling for Learning Model Predictive Control (L-MPC) with the initial condition Pk,o −1 = CΣ0C In most of bio-chemical batch processes there are no observations of the initial xk,0 0 Pk,0 0 = CΣ0C However, if these are available they can be specified. The Kalman Filter re- cursion is used for estimation of outputs or states when the measurement is notavailable on-line.
Benchmark Control Design
To implement the controller to the benchmark model a Learning Model Pre-dictive Control is used which was developed in the methodology. In thesubsection was presented the controller. In this subsection the isimplemented to the benchmark example (Penicillin Production). The controlobjective for this case study is batch reproducibility. The model used for thecontroller is the second best model estimated by in the previous sectionthe objective is to test the robustness of the controller instead the modelused is different from the plant.
The duration of the batch is 150 hours, the samples measurement (N ) for the outputs variables (yk,t) and inputs variables (uk,t) are each hour (N = 150). Theprocess starts with a batch phase where no feed is added, the initial substrateadded is sufficient to ensure biomass growth during this phase. When the sub-strate concentration is less than 5 g/L the feed starts. The initial conditions aregiven as k,0 = ⎣ 0 ⎦ + ⎣ where pk is a random variable sampled from a zero mean uniform distributionwith variance 1. The formulation of the control objective used is {Δuk,l,t}N−1 = arg min k,i }N −1 i=t+1 uk,i−1 − uk,i−2)Ti(uk,i−1 − uk,i−2) + ˆ xk,i t = Ai ˆ xk,i−1 t + BiΔuk,i−1 ek,i t = ˆek−1,i N − C ˆ uk,i−1 = Δuk,i−1 + uk−1,i−1 ymin ≤ ˆ ek,i t ≤ umaxi M.A. Alvarez, S.M. Stocks, and S.B. Jørgensen with uk,−1 = uk−1,0. Where u(1) is the feed flow rate and u(2) is the substrateconcentration in the feed. The output constraints for the process are ymaxi = ⎣ 10 ⎦ for the samples i = 1, . . , N . The measurements are: y(1): biomass concentration(g/L), y(2): substrate concentration (g/L), y(3): product concentration (g/L)and y(4): volume of the reactor (L). According to the control objective theweighting matrices were for this preliminary design selected as The control objective is batch product reproducibility using a previous batch as reference. For cultivation industries this type of control is interesting for re-ducing the influence of variation of initial conditions or substrate compositionon the batch product concentration. The above preliminary control design istested introducing variations on initial conditions (as specified in Eq. Fig-ure presents the open loop performance for simulation of 20 batches (Peni-cillin model), where the inputs vary within 15% of the nominal value. The mean Substrate feed concentration Fig. 6. Open-loop performance simulation of 20 batches. The inputs are perturbed
with an amplitude of 15% of the nominal value. Stochastic initial conditions.
Bioprocess Modelling for Learning Model Predictive Control (L-MPC) Fig. 7. Closed-loop performance of L-MPC control with stochastic initial conditions.
The control objective is to ensure product reproducibility.
Table 7. FIT for Pure Simulation Prediction of batches 6-7 for ARX model
Operation Mean Product Conc. Productivity Disturbance on initial conditions At 50 hours μmax decreases 30% product concentration of the 20 batches is 3.7 g/L. The batch with the highestproductivity has a final product concentration of 4.5 g/L. This batch is chosenas reference batch for the control design.
The closed-loop performance simulation for 20 batches, the learning MPC controller efficiently rejects the perturbation produced by variation on the initialconditions (see Figure Note especially the behavior towards the end of thebatches for product concentration. Table shows the mean value of the productconcentration without control and with control.
Another scenario for control testing is generating a disturbance during the batch, e.g. by decreasing the maximal specific growth rate (μmax) with 30%of the original value after 50 hours. This disturbance simulates a change inthe metabolism of the microorganism due to an inhibitory effect. Figure ) M.A. Alvarez, S.M. Stocks, and S.B. Jørgensen (a) Open-loop of the system with disturbance
(b) Closed-loop disturbance rejection
Fig. 8. (a):Open-loop with disturbance. The 30 % reduction of μmax is applied at time
50 hr. (b): Close-loop performance of a learning MPC control for disturbance rejection
intra-batch.
illustrates the open-loop performance with the disturbance in μmax at 50 hours.
Figure b) depicts the closed-loop performance of control, and showsthat even though the plant changes its behavior, the controller attempts to min-imize the error by reducing the feed rate while increasing the feed concentration.
Table shows that this is achieved quite well.
Bioprocess Modelling for Learning Model Predictive Control (L-MPC) Conclusions on Benchmark
On the benchmark simulation case it is demonstrated that it is possible to de-velop a data driven model from a limited number of batches using designed batchruns with perturbations in the input variables which are going to be used as ac-tuator variables in the control design. There are parameters in that can beused for model improvement such as, changing the number of batches for iden-tification, validation and testing (selecting the batches with input perturbationfor identification); the weighting matrix (λ) for Tikhonov Regularization; themodel type or maximum model order; the Model Purpose (a value between 0 to 1). The Benchmark also illustrates that it is possible todevelop a soft sensor for the enzyme activity, which otherwise normally is notavailable during the cultivation. Finally a preliminary control designis demonstrated to provide reasonable control both in case of intra-batch andinter-batch disturbances. As expected, even thought the plant changes its be-havior, the controller tries to reject the disturbance in spite of the model usedfor the controller is different from the plant. This allows to conclude that theseems to be a robust controller. Hence a suitable background has beenestablished for investigating the application of the methodology on anindustrial case.
Based upon the above benchmark investigation it was decided to investigateapplication of the methodology on a challenging fungal cultivation. Themain challenge stems from the fact that the theology of the cultivation is com-plicated by the growth of filaments wherein the main activity during growth andproduction actually resides.
Introduction to Pilot Plant
The Industrial Case Study is the production of an α-amylase called Fungamylusing the microorganism Aspergillus oryzae. The fermentation process is oper-ated following a standard recipe at the Novozymes Pilot Plant. The availableon-line and off-line variables are presented in Table The fermenter feed flow design has the physical connections illustrated in Fig- ure Two feed lines are available for the experimental cultivation to excite theinlet total feed flow rate and the substrate feed concentration to the reactor. AsFigure shows the inlet feed flow rate is the total feed (FT = Fd +Fw). The inletsubstrate concentration (Sf ) is total amount of substrate per mole of total feed.
The dosing tank is prepared before the main tank is inoculated according to the recipe; this tank is connected to the reactor to transfer the substrate, withthe flow rate Fd which is manipulated using control valve V2. The seed tankis filled with water after the inoculation has been transfer to the main tank;the water flow is manipulated with a control valve by a previously designedsignal. These two flows are connected before the fermenter, and the substrate is M.A. Alvarez, S.M. Stocks, and S.B. Jørgensen Table 8. Variables available from Pilot Plant
Dissolved O2 tension (DOT) O2 uptake rate (OUR) CO2 evolution rate (CER) Accumulated CO2 evolution Refractive index (RI) Respiratory quotient (RQ) Input variables
Feed dosing flow rate (Feedd), Measured L/h
Accumulated dosing Feedd flow
Dosing Feedd flow rate, Set Point Water Feed flow rate (Feedw), Measured L/hAccumulated water Feedw flow Feedw flow rate, Set Point Ammonia concentration Fig. 9. Flowsheet for Pilot Plant feed line connections
Bioprocess Modelling for Learning Model Predictive Control (L-MPC) Fig. 10. Connections between Pilot Plant computers
mixed with water and diluted to generate the perturbation in the inlet substrateconcentration and in the total flow rate.
The physical connections of the computer system are shown Figure All sensor signals are connected to the DeltaV System that communicated withthe DeltaV Server by an Ethernet. The PI Database Server stores all the dataacquired from the DeltaV Server using a "swinging door" algorithm. The hostPC directly collects data from the DeltaV Server every 10 seconds and this dataset is further prepared for analysis and model identification.
Given the above flow sheet (Figure it is important to carefully design of the perturbation signals in order to ensure that substrate concentration andtotal flow rate indeed are independently perturbed. Obviously this would not bethe case if such signals were applied directly to the valves V1 and V2. Thereforespecial consideration is given to the design of these two perturbation sequences.
Input Perturbation Design
The input variables FT and Sf are, as explained previously, closely related dueto the physical connections (see Figure . This fact is also evident from themass balances: T = Fd Sf = Sd where Sd is the Dosing Substrate. Therefore, the aim is to design Fd and Fw,which both are physically correlated, such that FT and Sf will be perturbednearly independently. The strategy is to linearize the equations and , M.A. Alvarez, S.M. Stocks, and S.B. Jørgensen and then decouple the equations. Defining the steady state values as ps =[F s (F s+F s )2 then the linearized model becomes: Define the covariance matrix for . such that their cross-covariance will The Lyapunov equation for the covariance of [FT Sf ] is (assuming no process P = ARAT where the A matrix is defined by equation −aF sd aFsw and σ21 and σ22 are the variance of Fw and Fd respectively. An R matrix is definedas: R = HHT The aim is to determine R such as the cross-covariance between FT and Sf is zero. Given the square root H matrix it is possible to determine ΔFd and ΔFw: where e1 and e2 are the perturbed signals. In this case these signals are uncorre-lated PRBS signal with period of 4 hours for e1 signal and 8 hours for e2 signal,and hence Fw and Fd become: Fw and Fd were designed following the above analysis. The reference signals forFT and Sf are ramp signals, both signals start at the beginning of the fed-batchphase. From these signals Fw and Fd are designed according to Eq. 58. Figure shows the resulting signals.
A correlation analysis for FT and Sf shows there is a cross-correlation of 0.0646, which indicates that both signals are almost independent. Thus with theabove linear based design it is possible to satisfactorily decouple the signals inspite of the physical coupling.
Bioprocess Modelling for Learning Model Predictive Control (L-MPC) Fig. 11. Top: Feed water flow rate signal; bottom: Feed dosing flow rate signal
For most microbial batch processes it is difficult to develop a model using firstprinciple modeling, thus, there is an interest to use black-box models. The aimin this case is to illustrate estimation and validation of a model using the methodology for an α-amylase cultivation of the Novozymes Pilot Plant. First,it is presented the data of the batches available for identification. Then, it isshown the estimation of the model using the validation of the model.
Data Selection: The eleven batches presented below are from production of
α-amylase in Aspergillus oryzae cultivations from the Novozymes Pilot Plant.
The available results in Table are produced in the same reactor.
Classifying the batch behaviors revealed that five of these batches fall out of the normal operation pattern of the process . From the available data set,only batches 1082, 1098, 1099, 1108, 1110 and 1111 are selected to be potentiallyuseful for model identification. After selection of the potentially useful batchesit is necessary to pre-process the data as explained next.
Variable Selection and Data Pretreatment for Modeling: It is important
to select variables for model estimation that contain sufficient information. The
M.A. Alvarez, S.M. Stocks, and S.B. Jørgensen Table 9. Variables available for each batch
Type 1082 1098 1099 1100 1101 1102 1103 1108 1109 1110 1111 Input perturbation (FT -Sf ) online NH3 flow rate Feed-back control Table 10. Carbon concentration in dosing flow. The concentration is constant during
each batch process, except for batches 1110 and 1111 where the carbon concentration
is perturbed during the batch as explained in section
Batch Calc. Carbon concent.
pH, back pressure and temperature are controlled, for that reason they are notused for modeling but their control actuator signals may be used. The air flowrate is kept constant during the cultivation and depends on the capacity of thecentral compressor. The remaining variables do contain information about theprocess and are used to estimate a model with the methodology.
The input variables to the process are: Feed flow rate and Substrate feed concentration. The decoupling design illustrated in Figure for the two inputvariables was used in batch 1111, while correlated perturbations were introducedin batch 1110. It is necessary to calculate the substrate concentration in the feed.
Table summaries the carbon concentration in the dosing for the batches usedfor model estimation. Figure depicts the total feed flow rate FT and substratefeed concentration Sf for all barches.
The data preparation is carried out for the variables than can be used for model estimation, first outliers is removed from the data. That is, data that isinconsistent with the normal operation by visual inspection is removed. At thispoint the dead time produced by delay response of sensors is removed. Then akernel smoother is used for re-sampling data to reduce the effect of noise and ofslow sampling interval for a few of the off-line sampled variables. Thereby highfrequency measurement noise is removed and missing data reconstructed. Caremust be taken not to remove system dynamics while filtering the measurementnoise. Not all the measurements are available for all the batches (see Table and depending of the variable different methods were used as explained next.
Figure shows the filtered variables for all the batches available for modelidentification.
Bioprocess Modelling for Learning Model Predictive Control (L-MPC) Substrate feed concentration (a) Total flow rate
(b) Substrate feed concentration
Fig. 12. Left: Total feed flow rate for all the batches; Right: Substrate concentration
in the feed for all the batches. Batches 1110 and 1111 were design with perturbations
in total flow rate and substrate feed. Only batch 1111 used the designed uncorrelated
perturbations.
The measurement of ammonia flow rate is not available for the batches 1082 and 1098. To estimate the ammonia flow rate, the accumulated ammonia that ismeasured for each batch was used. The accumulated N H3 is a measure of the to-tal amount of N H3 added to the fermentation tank and it is the integrated N H3flow rate. Therefore, the differentiation of the accumulated ammonia providesthe ammonia flow rate. This is relevant for batches 1082 1nd 1098.
The variables used for model estimation are: – Inputs: Total feed flow rate and Substrate feed concentration
– Outputs: N H3 flow rate, CER, DOT, RI, Viscosity and Enzyme activity
Model Estimation: Several different models have been estimated with
different combinations of outputs variables (the same inputs is using for estima-
tion of all the models [FT and Sf ]), different sets of batches and different values
of Model Purpose For the purpose of obtaining a model with reasonable
prediction ability to be used for the implementation of the control on the Pilot
Plant. The following mostly sequential methodology was applied:
1. Selection of output to be used for modeling.
2. Selection of batches used for estimation, validation and an independent set for testing.
3. Selection of model type.
4. Selection of maximum local model order.
5. Set the scalar for Model Purpose (ρ).
This section only presents the combination used to estimate the best model obtained for the Fungal Process. After many estimation attempts it was possibleto conclude that batch 1082 is not a good batch for identification. Table shows a comparison between two estimated models, one using batch 1082 and M.A. Alvarez, S.M. Stocks, and S.B. Jørgensen Fig. 13. Pilot plant data available for data driven model estimation
Bioprocess Modelling for Learning Model Predictive Control (L-MPC) Table 11. Comparison between two models, one with batch 1082 and the other without
this batch
Model Inputs Ouputs NH3 DOT CER Enz.
Nest:1110,1082,1099,1111 0.0756 0.1149 0.0868 0.1448 0.1089 0.0837 0.1012 0.074 0.1158 the other without this batch. It is interesting to note that without batch 1082, theestimated model (model1) improves its prediction. This indicates that batch 1082does not follow the real behavior of the process and its data does not providereliable information. Therefore, care should be taken when data is chosen formodel estimation. As only there are available six batches for model estimation,the batches 1110,1082,1099,1111 are used for estimation, batches 1098,1108 areused validation and testing. Figure shows the variables used for identificationincluding data from batch 1082.
According to the above results, it was decided to not include batch 1082 for model estimation. When the estimation of a model is using a black-boxidentification as the case of it is necessary to use data that containsa lot of information of the process. For that reason, a model was estimatedwith two more measured outputs (RI (Refractive Index) and Viscosity), as soonas these measured variables became available. Both new measurements provideinformation about biomass and product concentration during the cultivation.
Table shows the FITs of four models estimated.
Since only five batches are available for the estimation, these are divided into three batches for estimation and two for validation. Therefore, there is not anindependent set of batches for testing hence the batch set for validation is alsoused for testing. Calculation of the generalization error is, as for usedto continue the iteration until convergence where the Estimated GeneralizationErrors is less than the user specified value. This is supplemented with visualvalidation.
k − Δˆ θ∗ is the selected near optimal model θ∗ = ˆ V , W, Λ∗, n∗ The model model2 in Table which includes RI (Refractive index) is not as good as model1, while model3 which includes the viscosity measurement clearly M.A. Alvarez, S.M. Stocks, and S.B. Jørgensen Table 12. Model estimated with different combinations of outputs
Model Inputs Ouputs ARX Nest:1110,1111,1099 0.0837 0.1012 0.074 0.1158 Nest:1110,1111,1099 0.0896 0.1078 0.085 0.1382 0.107 Nest:1110,1111,1099 0.0754 0.0841 0.0788 0.0625 Nest:1110,1111,1099 0.0755 0.0873 0.0666 0.0873 0.0409 0.0838 0.0718 is better than the other two. However including both RI and the viscosity mea-surement visc2 improves the fit significantly, which indicates that these twovariables contain complementary information probably on substrate breakdownand on the influence of the microorganism filaments on the rheology.
Figure shows the validation of the best model visc2 on batch 1098. In the lower part of the figure the batch increment deviation between the model visc2and batches 1098 and 1108 clearly illustrates that even though the fit is quitereasonable for visc2 there are deviations in the estimated enzyme activity. Notehowever that there only are four measurements during a batch of this variable.
The other variables actually fit quite well, except during the shift from batch tofed-batch. The estimated model visc2 is being used in the following for KalmanFilter and control design.
To verify the influence of the perturbation on the substrate feed concentration using water to dilute the dosing two models were estimated. The models wereestimated using two inputs (FT and Sf ) and four outputs (N H3, DOT, CER andEnz.), the model structure used is with a maximum local model order oftwo and with a value of = 0.618. The first model was estimated using batches1099 and 1098 for identification and batches 1108, 1110 and 1111 for validation,for testing of the model was used the same batch set as for validation. Table shows the results of the fits for each model. The second model estimated wasused batches 1110 and 1111 for identification and batches 1099, 1098 and 1108for validation and testing. Figure presents the validation of both models. Themost pronounced effect is seen on the DOT (and the ammonia) fit which is poorwhen there are no deliberate perturbations (model test1) on the input variables.
Also the overall fit is significantly better with perturbations for (model test2)which has perturbations on the inlet flow variables.
Bioprocess Modelling for Learning Model Predictive Control (L-MPC) Pure Simulation Prediction (a) Batch 1098
Batch incremente: Batch4 − Batch5 Pure Simulation Prediction (b) Batch increment: Batch 1098 - 1108
Fig. 14. (a): Validation of the model visc2 using batch 1098. (b): Batch increment
between batches 1098 and 1108.
M.A. Alvarez, S.M. Stocks, and S.B. Jørgensen Pure Simulation Prediction (a) Batch 1108 - model test1
Pure Simulation Prediction (b) Batch 1108 - model test2
Fig. 15. (a): Validation of the model test1 to simulate batch 1108 using as previous
batch 1098. (b): Validation of the model test2 to simulate batch 1108 using as previous
batch 1110.
Bioprocess Modelling for Learning Model Predictive Control (L-MPC) Table 13. Model estimated with different combinations of outputs
Model Inputs Ouputs ARX 0.195 0.3407 0.1124 0.1322 0.1947 0.0806 0.1464 0.1168 0.1429 0.1367 Fig. 16. Closed-loop Scheme for the Fungal Cultivation. It is identifying the controller
block, process block and Kalman Filter block conections.
Having developed a model with reasonable predictive power, then the follow- ing section details the implementation of a soft sensor for enzyme activity on asimulation and of the design and closed-loop simulation of control of theFungal Cultivation (Figure .
Pilot Plant Control Simulation
The model obtained (visc2) is intended for the designing a controllerfor implementation on the Pilot Plant. However first the designed controller istested in simulation where the same estimated model, as used for control design,is used to simulate the process. However since the enzyme activity is an off-line (seldom) measured variable it is necessary to estimate it. For that purpose aKalman Filter is designed as a soft sensor to estimate the values of these outputs.
The basis for the filter design is also the estimated model. Below the open-loopmodel is simulated to investigate the tuning of the Kalman Filter. Once theKalman Filter is working, the closed loop process is simulated to investigate thecontroller tuning.
M.A. Alvarez, S.M. Stocks, and S.B. Jørgensen The control objective is batch reproducibility and the control design method- ology is as specified in section The following formulation is used: {Δuk,l,t}N−1 = arg min k,i }N −1 i=t+1 uk,i−1 − uk,i−2) ˘ Ti(uk,i−1 − uk,i−2) + ˆe xk,i t = Ai ˆ xk,i−1 t + BiΔuk,i−1 ek,i t = ˆek−1,i N − C ˆ uk,i−1 = Δuk,i−1 + uk−1,i−1 ymin ≤ ˆ ek,i t ≤ umaxi Soft sensor Design: During the Fungal Cultivation the measurement of en-
zyme activity is not available on-line, as this variable is an off-line measure-
ment where it takes days to obtain the biochemical analysis result. Therefore,
a Kalman Filter is used to estimate the enzyme activity. The Kalman Filter is
designed using the estimated fermentation model. The process variable measure-
ments are assumed containing a normally distributed measurement noise with
mean 0 and variance 0.1. The Kalman Filter tuning used is: The diagonal ele-
ments of P (Sate Covariance Matrix) are one, the diagonal elements of R (Noise
Covariance) are 0.01. The Disturbance Covariance matrix is estimated during
model estimation by the Software. Note that this disturbance covariance
matrix estimate is most useful for obtaining a suitable tuning of the Filter. Fig-
ure illustrates the performance of the Kalman Filter on batch 1098 with noise
on the output variables. Note that even though the enzyme activity is an off-
line measurement, it is indeed possible to estimate the variable on-line using the
Kalman Filter.
Closed-loop Simulation: A closed-loop simulation of the Fungal cultivation is
carried out using the model estimated and a control design as specified
below and detailed in section The control objective is batch reproducibility
and the batch used as control reference is batch 1111. This batch was chosen
because the DOT is kept high with a low viscosity, this enables increasing the
inputs without the problem of lack of oxygen, and furthermore, the enzyme
activity is higher. The output constraints for the process are
i = ⎢ 15 ⎥ Bioprocess Modelling for Learning Model Predictive Control (L-MPC) S T concentration Fig. 17. Green: Process Measurement with noise; black: The Kalman Filter estimated
output; red: The real output measurement of batch 1098.
for i = 1, . . , N and y(1) is N H3 flow rate (g/h), y(2) is DOT (% of saturation),y(3) is CER (mole/h) y(4) RI, y(5) is the viscosity (cP) and y(6) is the enzymeactivity (FAU/g). Following the control objective the weighting matrices weredesigned as The results are shown in Figure estimates the model for the dif- ference between two batches, as is yk = Δˆ yk + yk−1. Therefore, to simulate the plant using a model, it is necessary to use a previous batch, in this casesbatch 1098 was used. The control objective is product reproducibility. Duringthe batch phase the controller is open since no feed is added, the criterion tostart the feed is when the pH is above 6.3. For the simulation, there is no mea-surement of pH hence time is used to start the feed because all the batches startmore and less at the same time. Thus at time 35 the control loop is closed andfeed starts. Figure shows that the enzyme activity is tracked rather well, eventhough the fed batch is started a little later than for the reference batch. How-ever this promising controller needs a more careful investigation of the tuningbefore being implemented in practice.
M.A. Alvarez, S.M. Stocks, and S.B. Jørgensen Fig. 18. Closed-loop simulation. The plant is simulated using the model visc2 with
batch 1098 as previous batch while batch 1111 is the reference batch for the control
trajectory.
Discussion and Conclusions
Computational intelligence is a branch of science dealing with problems thatcannot be solved using only effective computational algorithms. Computationalintelligence combines elements of learning, adaptation, and modelling to developsolutions that are, in some sense, intelligent. The modelling methodology pre-sented in this chapter may be viewed as intelligent control ) in the sense ofqualitative control where a qualitative model is established based upon operationexperience in the form of knowledge about the batch plant purpose combinedwith modelling the behaviour between samples. The plant purpose knowledge isessential for for structuring the more detailed data driven modelling of the be-haviour between samples. The type of models selected for modelling the behaviorbetween sample points could be, e.g. based on neural nets. However It is impor-tant to stress that even though the behaviour between measurement points isapproximated with linear models then the behaviour is indeed timevarying sincethe parameters may vary between measurement points. Thus the selected linearmodels seem justified provided the measurement points are sufficiently close.
Hence the methodology presented and applied in this chapter, is intelligent inthe sense that intelligent modeling, which could have used fuzzy logic or neu-ral networks instead of model structure, is applied to model thebatch behavior based upon our knowledge of batch structure.
Bioprocess Modelling for Learning Model Predictive Control (L-MPC) Due to the repetitive nature and the finite horizon of the process, the mea- surements are collected from a grid of sample points in time. This means thata sample of all possible measured variables is collected as often as possible, i.e.
at each sampling period. With these sampled data, it is possible to model theevolution between two consecutive sample points in a batch with local mod-els, which are labeled grid-point models. This modeling is made possible since ateach sample time there are several measurements available due to the repetitionof batches. The Methodology estimates a model which contains alarge set of smaller models describing the process behaviour from batch tobatch, while the model describing the behaviour within the batch is The methodological choice in this paper is dictated by our knowledge of three functional requirements or facts. The first of these requirements is the desire to beable to reject inter and intra-batch production process disturbances. The secondis the fact of limited available knowledge of the microbial physiology wherethe regulatory networks still are largely unknown. Consequently a data drivenmodeling methodology has to be selected. The third requirement originates fromconsidering the control methodology where it is considered mandatory to beable to execute control actions fairly often during cultivation to prevent thecultivation from entering into operating regions where irrecoverable microbialstress responses may be triggered which might lead to opening of new pathwayswhere substrate may be lost to other purposes than production of the desiredproduct. This third overriding consideration lead to the choice of a discrete timemodel. The discrete time model is developed using the approach from acombination of existing operating data and designed fed batch runs with inputperturbations in order to obtain informative data for the modeling for control.
Furthermore the importance of using measurement variables that contain mostrelevant information about the process is demonstrated.
The modeling results of this paper demonstrates that it is actually feasible to obtain a reasonable cultivation model from a few designed production exper-iments. In fact it turned out when the experienced personnel saw the responsesfrom the input perturbations they were enthusiastic and wanted to employ thesefor other purposes as well.
The obtained discrete time model is demonstrated to be directly applicable for providing a soft sensor signal of the measured variables, such as enzyme activityin this particular case. The tuning of this filter is greatly facilitated by utilizingcovariance estimates from the estimation.
The control contribution presents how asymptotic convergence of Iterative Learning Control may be combined with the closed-loop performance of ModelPredictive Control to form an asymptotically stable optimal controller labeledLearning Model Predictive Control for ensuring reliable and reproducible opera-tion of batch and periodic processes. Consequently this controller type may alsobe used for optimizing control, . The data driven approach obviously restrictsthe usability of such a model to similar cultivations of the same microorganism,however if the product is a bulk product which is produced relatively often, thenmany batches will be able to benefit from a soft sensor and control design, such M.A. Alvarez, S.M. Stocks, and S.B. Jørgensen that there is a sound economical basis for such an investment. The benefit isobviously even larger if a high value product is produced in large amount.
The Methodology provides a methodology for black-box model estima-tion complemented with the for handling control both within the batchand between batches. The methodology has in this chapter been ap-plied to an Aspergillus oryzae cultivation, where the produced enzyme is usedto speed up the degradation of starch to glucose. Since the microorganism isa filamentous fungus this production process is very challenging. However it isdemonstrated that the GoLM methodology can be used for model estimationwhen only few batches are available. However, it is important that the batcheshave been designed for identification, this means, that the data set recorded hasrelevant information on actuator movement. The importance of the informationcontent is demonstrated on the benchmark case study through applying differentperturbation amplitude on the inputs.
With the obtained model a was developed and implemented. The controller was tested under two different scenarios; the first scenario simulatesdisturbances in the initial condition, which is a typical disturbance in BatchFermentations. The control objective is batch reproducibility, the control rejectsthis disturbances and a model obtained with the same productivity as specifiedby the reference batch. Such a reproducibility is a significant advantage whenthe cost of loss for the variations on the final product is significant. In the secondscenario an intra-batch disturbance occurred. In this case the robustness of thecontroller was tested during the simulation where the maximal growth rate waschanged, this. The controller effectively minimized the error, even though themodel used for the controller did not represent the plant due for the parameterchange.
The experimental design for the industrial cultivation illustrates the impor- tance of perturbing the process input. An amplitude of 25% of the nominal valueyielded informative measurement data, as realized during two cultivations. Con-catenating these batches to previous unperturbed batches it was possible to havesufficient data for model identification. From the previous batches it wasmost important to eliminate batches with a qualitatively different behaviourmainly using classification.
Based upon the discrete time model the development of a soft sensor for the measured variables including enzyme activity is demonstrated. This sensorrepresents perhaps the first direct application of the model. The model mayalso be used to optimize the Batch Operations Model for the cultivation. Andfurthermore the model is directly applicable for design of a controllerfor handling intra- and inter-batch disturbances as illustrated in simulation. Itis important to emphasize that even though the model of the process used forthe control may not be perfect, then the controller can iteratively minimize thiserror to reject intra and inter batch disturbances and ensure stable learning bysuitable design of the control weighting.
Bioprocess Modelling for Learning Model Predictive Control (L-MPC) ARMAX AutoRegressive Moving Average models with eXogenous
ARX AutoRegressive models with eXogenous
Gρ
Estimated Generalization Errors GoLM Grid of Linear Model
Iterative Learning Control L-MPC Learning Model Predictive Control
Linear Time-Invariant LTV Linear Time-Varying
MPC Model Predictive Control
NCO Necessary Conditions of Optimality
OSA One Step Ahead
PRBS Pseudo Random Binary Signal
Tikhonov Regularization wLS weighted Least Squares
1. Bajpai, R.K., Reuss, M.: A mechanistic model for penicillin production. J. Chem.
Technol. Biotechnol. 30, 332–344 (1980) 2. Bajpai, R.K., Reuss, M.: Evaluation of Feeding Strategies in Carbon-Regulated Secondary Metabolite Production through Mathematical Modeling. Biotechnologyand Bioengineering 23, 717–738 (1981) e, D.: Optimal and Reproducible Operation of Batch Processes. Technical University of Denmark, Department of Chemical Engineering. Ph. D. Thesis (2005) 4. Flores-Cerrillo, J., MacGregor, J.F.: Control of Batch Product Quality by Trajec- tory Manipulation Using Latent Variable Models. Journal of Process Control 14,539–553 (2004) 5. Kristensen, N.R.: Fed-Batch Process Modelling for State Estimation and Optimal Control. Technical University of Denmark, Department of Chemical Engineering(2003) 6. Lind, M.: Status and Challenges of Intelligent Plant Control. A. Rev. Control 20, 7. Ljung, L.: System Identification: Theory for the User. Prentice-Hall, Englewood 8. Petersen, N.: Multivariable Modeling for Control of Industrial Fed Batch Culti- vations. Department of Chemical ´ Engineering, Technical University of Denmark, M.Sc. Thesis (2006) 9. Rawlings, J.B., Bonn´ e, D., Jørgensen, J.B., Venkat, A.N., Jørgensen, S.B.: Un- reachable Setpoints in Model Preditive Control. Transactions of Automatic Control(2008) M.A. Alvarez, S.M. Stocks, and S.B. Jørgensen 10. Srinivasan, B., Bonvin, D.: Dynamic Optimization under Uncertainty Via NCO Tracking: A Solution Model Approach. In: BatchPro - Symposium on KnowledgeDriven Batch Processes, Kipparisidis, C, pp. 17–35 (2004) 11. Srinivasan, B., Biegler, L., Bonvin, D.: Tracking the Necessary Conditions of Opti- mality with Changing Set of Active Constraints Using a Barrier-Penalty Function.
Comput. Chem. Engng. 32, 572–579 (2008) Akesson, M., Hagander, P., Axelsson, J.P.: Probing Control of Fed-Batch Cultures:Analysis and Tuning. Control engng. practice 9, 709–723 (2001)

Source: http://www.fiec.espol.edu.ec/publicidades/2010/cd/DATA/papers/Bioprocess%20Modelling%20for%20Learning%20Model%20Predictive%20Control%20(L-MPC).pdf