- Email: [email protected]
Development of learning control for reproducible and high quality operation of batch processes
IFAC
Copyright © IFAC Dynamics and Control of Process Systems, Jejudo Island, Korea, 200 I
c:
0
[>
Pu blications www.elsevier.comllocate/ifac
DEVELOPMENT OF LEARNING CONTROL FOR REPRODUCIBLE AND HIGH QUALITY OPERATION OF BATCH PROCESSES Dennis Bonne * , 1 Sten Bay
*
J~rgensen * ,2
Department of Chemical Engineering, Technical University of Denmark, DK-2BOO Lyngby, Denmark
Abstract: In this contribution, a framework for modeling and control of batch processes has been extended with the effect of initial conditions on batch evolution and with the inclusion of measurement noise. The imbedded model state estimation algorithm has been stabilized and extended with the ability to process noisy measurements. A scheme has been developed for identification of high dimensional models and it has been shown that models with immense dimensions can be obtained from sparse, noisy data. The extended framework and identification scheme were validated through simulated trajectory tracking of fed-batch fermentations . Copyright © 2001 [FAC
Keywords : Batch process, modeling, identification, model predictive control, iterative learning control, yeast fermentation
1. INTRODUCTION
dition , implementation of control will not only increase reliability, but also enable process optimization. The modeling and control framework presented in this paper is designed to reject the two first mentioned types of variation.
With an increasing industrial demand for flexible and specialized production methods, batch and semi-batch processes are becoming ever more popular. Control of batch operation, nonetheless, has received little attention in academia. Most batch processes are operated according to off-line specified operation plans, which typically specify reference trajectories for physical and chemical properties. Batch and semi-batch processes however, often deviate from normal/ desired operation due to three types of variation: Variation that occurs during a batch run, variation between batch runs, and variation that reappears after a certain period. These variations from normal operation cause operation scheduling and plant design to be conservative. Consequently, to reduce the conservatism of present operation scheduling and plant design , feedback control is required to reduce variations from normal operation. In ad1 Corresponding author; fax: +45 4593 2906 e-mail: db@kt .dtu .dk 2 fax: +45 4593 2906 and e-mail : [email protected] .dk
2. MODEL DEVELOPMENT The following batch modeling approach is based on a framework (Lee et al., 2000) where nonlinear and time-varying dynamic behavior is approximated with a set of simple, Linear Time-Invariant (LTI) models. Lee et al. (2000) suggested using a set of Finite Impulse Response (FIR) models and Gregersen (2001) suggested using a set of AutoRegresive models with eXogenous (ARX) inputs. Here the model may be set up using either FIR models or ARX models . Operation modes of batch processes will vary with time. Therefore, the number of measurements (which require trajectory tracking) and actuators may vary with time . Often it will be necessary to include additional measurements for monitoring and estimation purposes, and to measure the quality of the
and
383
the resulting combined model is:
end-product. Consequently, it seems reasonable to define the following variables:
ek+1 = ek -
• Control variable, u(t) E Rnu(t). • Controlled output variable, y(t) E Rny(t), which requires tracking control. • Secondary output variable, s(t) E Rns(t), used only for monitoring. • Quality variable, q ERn., only measured after a batch.
ek = ek
y:=[yT(l)
]T , YT (2) ... yT(N)
T
d:= [d (l)
o
ek+1
(1)
-
GY ~Uk+1 -
e~:=y-y=e~+v~,
One major drawback of the proposed modeling approach is the immense dimensionality of the resulting set of models. E.g., assuming causality a process with one output and one input sampled 100 times during a batch will be parameterized by 5050 parameters using FIR models and using ARX models 198 model orders and 594 parameters assuming all model orders set equal to 3. In practice, this immense dimensionality will render any identification problem ill-conditioned. Fortunately, there is some structural information in the resulting set of models G. Firstly, by causality the upper triangular is 0, secondly, the model parameters time variation can reasonably be assumed smooth (L 1 ), i.e., the diagonal elements in G will vary in a smooth manner, and thirdly, the impulse responses can reasonably be assumed smooth (L 2 ), i.e. , the column elements in G will vary smoothly. Regularization utilizes such a priori knowledge, reducing the amount of data needed, by introducing constraints or penalties in order to attract excessive degrees of freedom toward reasonable values. Unfortunately, regularization will result in biased estimates, i.e., there will be a tradeoff between variance and bias. The optimal tradeoff is found by selecting the regularization weights Ai that minimize an average prediction error over the validation batches. The regularized system:
+ Mni(Yini).
(3)
Grni~Yini,k+l
+ w~, (4)
Uk+1 - Uk and ~Yini,k+1 = The matrices GY and Grni are linear system approximations and e~ is the part of e~ that reappears in e~+l' Both w~ and v~ are assumed zero-mean i.i.d. (independent and identically distributed) sequences with respect to k.
where ~Uk+l
Yini,k+1 -
eZ + ~~,
3. IDENTIFICATION
The nonlinear model (3) may be approximated by a set of linear models. In the present paper N(u , d) will be approximated with one linear model for each sample instant. Let y be the specified and preferably optimal output reference trajectory and let k be the batch index. Then the output error trajectory, e Y , is defined as:
= e~
(7)
(2)
:= y(O) .
Note, not all initial conditions are measurable and/or physically meaningful. The output sequence is in general related to the input sequence, initial condition, and disturbance sequence by a nonlinear model:
e~+l
(6)
where vg is assumed a batch wise white disturbance term , z~ is the measurement of e~ , and ~~ is assumed batch wise white measurement noise.
where N is the batch length. Let the initial condition be defined as:
y = M(u, Yini , d) = N(u, d)
= ek0 + Vk0
z~ =
]T , dT (2) . .. dT(N)
Yini
+ Wk,
To enable usage of a Kalman filter for estimation of the initial condition eO, the initial condition is modelled with a simple batch wise random walk model:
The input, output and disturbance sequences over the entire batch are defined as: u T (l) ... u T (N-1) ]T ,
- G ini eZ+ 1
2.3 Model for Initial Condition
2.1 Model for Controlled Output
u:= [uT(O)
G~Uk+1
+ Vk·
=
Yini,k·
2.2 Combined Model
b reg = AregCf',
By applying the same modeling approach to the secondary and quality outputs and defining:
(8)
is solved in a least squares (LS) sense. The vector Cf' is the unknown model parameters, while Areg contains the system inputs and initial conditions (and outputs for ARX models) and breg contains
(5)
384
[!:::] = [~ ~] [!:::=~] - [~:=~] ~Uk(t Zk(t) =
y(/)
[0
Ht] [ek,t]
ek,t
- 1),
+ ~(t),
(11)
where t = 1, ... , N, Zk(t) is the measurement of the output error ek(t) at time t, and ~(t) is assumed zero-mean i.i.d. measurement noise at time t. Note, this state-space formulation does not seem observable, i.e., e.k,t is not observed in Zk(t) . However, remembering that ek = e.k+Vk, e.k ,t can in fact be observed through ek,t with a filter. At the start of a batch the states in (11) satisfy:
1(0)
Fig. 1. Block diagram of an Iterative Learning batch Control (ILC) algorithm in a Model Predictive Control (MPC) framework with combined output trajectory tracking and end-product quality control. the system outputs. The regularized system (8) is given as:
[ o~] [>.11] =
c{>,
(9)
>'2L2
Note that the structure of Ht and ek(t) depend on which measurements are available at sample t . Let €k,t+mlt be the optimal estimate of the error trajectory of the kth batch given measurements up to time t and m future control moves. Then the estimate €k ,t+mlt is given by:
where the exact system set up from Nk batch data sets:
(10)
€k,t+mlt = €k ,tlt - G'(' ~uk.t·
In literature, several approaches to selection of scalar optimal regularization weights have been presented. However, seemingly no work has been done with vectorial regularization weights based on model properties. For further discussion of how to solve the regularized LS problem (8) the reader is referred to Hansen (1996). Furthermore, usage of high quality data for the identification will decrease estimation variance. The quality of batch identification data can be determined by statistical comparison to a normal reference model obtained through multi-way analysis. For further discussion of high quality or normal batch identification data the reader is referred to Gregersen and J0rgensen (1999) and Nomikos (1996).
(13)
The optimal input movements are determined by employing the standard quadratic objective function known from conventional MPC control.
where Qt and R t are weighting matrices. The analytical solution to (14) is given by:
~uk.t = [G'('T QtG'(' + Rtt1G'('T Qt€k,tlt , (15) where only ~Uk(t) from ~uk\ is implemented. In practice input signals and their rate of change must be constrained to ensure non-provoking or system friendly and physically meaningful inputs. If possible the inputs should be constrained such that the process avoids regions of instability increasing the level of safety in the event that the controller should set out. Furthermore, the predicted outputs may need to be constrained to meet additional safety requirements and ensure operation in a physically meaningful operation region. Therefore, (14) is extended to:
4. MODEL-BASED ITERATIVE LEARNING
CONTROL Based on the newly obtained batch transition model, Iterative Learning batch Control (ILC) in a Model Predictive Control (MP C) framework can be applied for combined trajectory tracking and end-product quality control (Lee et al., 2000). Model predictive control has been comprehensively covered in literature, e.g. by Mayne et al. (2000). The block diagram of the control algorithm with the extensions introduced in the present paper (i.e., initial conditions and measurement noise) is shown in figure 1. From eq. (6) the following periodical, Linear Time-Varying (LTV) state-space model is obtained for control design:
Win
~Uk , t,Ek , t
[-2 {~uk.r Qt~uk.t + 2R; ~uk.t 1
+ELStEk,d], S.t.
At~uk.t::; €k ,t ~
385
0,
Bt(Ek,t),
(16)
'R;
where Qt = Gr;'T QtGr;' + R t and = -e.LltQtGr;' . By proper configuration of A and B both the outputs and the inputs may be constrained. The constraint softening slack variable E ensures feasible solutions and St is the corresponding weighting matrix. Since the secondary output is not controlled, the weighting matrix Qt has the following shape: Qt
=
Qi 0 [ o
1
0 0 0 0 0
Ql
Gt -
Pk,t1t-l] = [ ~k,tlt-l P k,tlt-1 P k,tlt-1
]
(20)
.m
T T
+ GiniPeoklkGini + R w ,
= Pk-1 ,NIN +GiniPeoklkG'fni + Rw
+ Rv .
In the above P is the prediction error covariance of the error trajectory, P is the prediction error covariance of the part of error trajectory that reappears, P is the cross-covariance between e and e, 'P is the covariance of the innovations, K is the Kalman gain, lC is the normalized Kalman gain, and Rn,t is the measurement noise covariance at time t . Note that the noise sequences are assumed both mutually independent and independent of the error sequences at time t . For further discussion of numerically stable implementations of Kalman filters the reader is referred to Grewal and Andrews (1993) and Kailath et al. (2000) . 4.1.1. Estimation of Initial Conditions The estimated initial condition eZ,k is obtained from: AO
ekl k =
(1
0 - "'k ) eAOk- 1Ik-1 + "'kZk>
+ Rvo) x (PeOk-1Ik-1 + Rv o + Rn o) -1, (1 - "'k) (PeOk-1Ik-1 + Rvo) ,
"'k = (PeOk-1Ik-1
PeO k1k =
(21)
where PeO k1k is the prediction error covariance of the initial condition, "'k is the Kalman gain, Rvo is the covariance of v O, and Rno is the measurement noise covariance.
1
P k,t-1 It-l]
' 5. SIMULATION RESULTS
~k,tlt] = [~k,tlt-1 ]
To investigate the applicability of the proposed modeling approach and controller design in biochemical industry a yeast producing fed-batch fermentation process was studied. To simulate fed-batch fermentation of yeast a biochemically structured model presented by Lei et al. (2001) was used. The model describes aerobic growth of Saccharomyces ceravisiae on glucose and ethanol
ek,tlt-l
-K k,dHte.k ,tlt-l - zdt)], 1t-1] lC k,t -- [Pk,t p tHT ' k,tlt-1 'Pk ,t = H t P k ,tlt-1 H ; + Rn ,t,
[gini] e2 ,k'
-
Pk,OIO
and the updating part is,
[ ek,tlt
+ Rn,t] -1 ,
= Pk-1 ,NIN + GiniPeoklkGini + Rw,
A
(18)
P k,t-1It-1 P k,t-1It-1
T
x [H t P k ,tlt-1 H ;
Pk,OIO = Pk-1 ,NIN
t:.Uk(t - 1),
[~k,t-1It-l
-
+ ICk,tKk,t'
K k,t = [P k,t 1t-1 ] H; P k,tlt-1
-
The optimal (i.e., minimum variance) estimates of the model states are obtained from a Kalman filter. Here the Kalman filter can be thought of as being made up of a set of Kalman filters or as an Extended Kalman Filter (EKF) without subsampling. By defining Rw and Rv as the covariance matrices of Wk and Vk, respectively, the predictive part of the stabilized Kalman filter (Thornton and Bierman, 1978) is,
1
-T
+Kk,tICk ,t
Pk,OIO
4.1 Kalman Filter
_ [G t -
Pk ,t1t-1] P k,tlt-1 P k,tlt-1
~k'OIO] = [~k-1'NIN] _ [ ek,OIO ek-l,NIN
The presented control technique transfers the refined control error sequence as initial condition for the next batch. Furthermore, it computes the input sequence as a change from the previous batch based on the control error prediction. This use of information from previous batches provides the control algorithm with an offset elimination action along the batch index. Thus, through this action the control algorithm is expected to eliminate offset errors both from disturbances that reappears in subsequent batches, and from model bias, i.e., the disturbances represented by w. A convergence prove for a batch MPC controller based on statespace model (11) can be found in previous work by Lee et al. (2000). This proof also holds for this extended algorithm.
[!:: :::=~] = [!:::=:::=~]
[~k,tlt-l
with the initial conditions,
(17)
.
~k,tlt Pk,t1t] = [ Pk,tlt Pk ,tlt
(19)
lCk ,t = 0.5K k,t 'Pk,t - lCk,t,
386
• The model validation was based on three normal batch data sets.
and focuses on the pyruvate and acetaldehyde branch points, where overflow metabolism occurs when the growth changes from oxidative to oxidoreductive. The fermentation was fed with substrate and the feed flow rate was manipulated as the control variable. In the simulation the ethanol concentration (EtOH), Oxygen Uptake Rate (OUR), and Carbon dioxide Evolution Rate (CER) were measured on-line, and it was assumed that there was no time delay of the measurements. From the on-line measurements OUR and CER the Respiratory Quotient (RQ) was calculated: RQ = CERjOUR. To characterize this fedbatch fermentation the following dimensionless measures were defined:
To facilitate the identification task a highly exciting input sequence was chosen. The input sequence was constructed by adding a Pseudo Random Binary Sequence (PRBS) to the nominal input trajectory. At each sample the PRBS was uniformly scaled by maximum ±1O% of the nominal trajectory. Note, highly exiting input sequences are not system friendly and thus seldom applicable in practice. For further discussion of input sequences see Pearson (1998). Since the basic assumption of this inter batch modeling approach is that trajectories vary only slightly from batch to batch, the validation batch data sets were generated with input trajectories as those used for identification, i.e., data sets similar to the identification data, but not used in the identification. From table 1 it can be seen that quite reasonable models can be obtained from noise free data for all the outputs with relative prediction fits of approximately 2 x 10- 2 . Furthermore, as indicated in table 2 models of similar quality can be obtained from noisy data (the noise level is based on laboratory experiments and industrial data) without pre-filtering. One of the three model cross-validations from the noisy identification is shown in figure 2.
• Quality == final yeast concentration minus final ethanol concentration x 1000 and substrate concentration x 200 normalized by unit concentration (1 9/l) . • Yield == produced yeast divided by fed substrate on a mass basis. The weights on ethanol and substrate concentration in the quality measure could be based on costs of downstream processing relative to the sales price of yeast, but here they are set such that variations in yeast, ethanol, and substrate concentrations are equally represented in the variation of quality. 5.1 Identification
5.2 Control
For simplicity the fermentation was modeled with EtOH, OUR, CER, and RQ as output variables and no quality variables. RQ was treated as a measurement to test if the dynamic behavior could be modelled with the proposed mode ling approach. To study the applicability of the proposed modeling approach to processes where initial conditions are not obtainable, the dynamics in G ini were discarded. A set of FIR models was identified in the following scenario:
The control simulations were designed to accentuate the iterative learning batch control feature and study the relation between trajectory control and reproducibility. The performance of a Single Input Multiple Output (SIMO) ILC algorithm without quality control was validated by randomly perturbing the initial conditions with max ±1O% in each batch run and introducing a persistent bias in the input trajectory by lowering the feed substrate concentration with 10%. The benchmark which the performance was tested
• The measurement noise was approximated by white noise with a standard deviation equal to one third of the maximum noise level. Maximum measurement noise: EtOH, ±3%; OUR,±2%;CER,±2%;RQ,±2%;u,±1%. • The initial conditions were uniformly perturbed by maximum ±1O%. • The identification was based on ten normal batch data sets. • The identification was optimized by testing possible combinations of the regularization weights ),1 and ),2 and choosing the combination that minimized the mean prediction error when cross-validated. The optimal regularization weights are shown in tables 1-2 along with their respective mean prediction fits.
Table 1. Optimal regularization weights and relative prediction fit for noise free identification. Output variable ~EtOH ~OUR
~CER ~RQ
'\1
'\2
fit
1.18 0.05 0.05 0.02
0.04 0.01 0.01 0.01
1.5746x 10 =-2" 1.9223x1O ·2 1.7818x1O =-2" 1.6298 x 10-"""7"
Table 2. Optimal regularization weights and relative prediction fit for noisy identification. Output variable ~EtOH ~OUR
~CER ~RQ
387
'\1
'\2
fit
1.18 0.05 0.05 0.26
0.05 0.01 0.01 0.02
1.5708 x 10 -2 2.4224x1O= 2.0479x 10 -2 2.3902x1O =2
batch evolution and with the inclusion of measurement noise. The modeling framework results in models with immense dimensions and hence an identification scheme has been developed utilizing structural information through regularization. Applying this identification scheme, high dimensional models can be obtained from sparse, noisy data. The model state estimation algorithm included in original framework has been stabilized and extended with the ability to process noisy measurements. The extended modeling and control framework and identification scheme were validated through simulated trajectory tracking of fed-batch fermentations with satisfactory results. Hence, the present extended framework for modeling and control of batch processes exhibits significant potential for industrial implementation.
~~] i]:::~ !j;:;~ ..,
~
o
I
2
3
..
5
•
7
e
,
10
o
I
2
3
..
5
11
7
IS
•
10
0
1
2
3
..
5
,
7
•
•
10
OO:f ~·050~---7-,--!-,-'----!3--.C-----!-'--!-,---!,:---,C---'-,---.J,O Time (h)
Fig. 2. Validation of noisy identification in batch incremented variables. True batch difference, solid line; predicted batch difference, dotted line.
i~b: : Z\ ::J !~t::52:1 !::~ ,
,
3
..
5
6
7
8
10
1
2
.,
..
5
•
7
a
10
1
2
3
•
5
6
7
a
7. REFERENCES
Gregersen, Lars (2001) . Monitoring and Fault Diagnosis of Fermentation Processes. PhD thesis . Technical University of Denmark. Gregersen, Lars and Sten Bay J0rgensen (1999). Supervision of fed-batch fermentations . Chemical Engineering 10umal75, 69- 76. Grewal, Mohinder S. and Angus P. Andrews (1993) . Kalman Filtering - Theory and Practice. Prentice Hall. N. J . Englewood Cliffs. Hansen, Per Christian (1996). Rank-Deficient and Discrete Ill-Posed Problems. Polyteknisk Forlag. Lyngby. Kailath, Thomas, Ali H. Sayed and Babak Hassibi (2000). Linear Estimation. Prentice Hall. N. J . Englewood Cliffs. Lee, Jay H., Kwang S. Lee and Won C. Kim (2000) . Model-based iterative learning control with a quadratic criterion for time-varying linear systems. Automatica 36(5), 641-657. Lei, F ., M. Rotb011 and S. B. J0rgensen (2001) . A Biochemically Structured Model for Saccharomyces cerrevisiae. Accepted for Journal of Biotechnology. Mayne, D. Q., J. B. Rawlings, C. V. Rao and P . O. M. Scokaert (2000). Constrained model predictive control: Stability and optimality. Automatica 36, 789--814. Nomikos, Paul (1996) . Detection and diagnosis of abnormal batch operations based on multi-way principal component analysis. ISA Transactions 35(3), 259-266. Pearson, Ronald K. (1998) . Input sequences for non linear modeling. Nonlinear Model Based Process Control. Proceedings of the NATO Advanced Study Institute pp. 599-62l. Thornton, C. L. and G. J . Bierman (1978) . Filtering and error analysis via the U DUT covariance factorization . IEEE Transactions on Automatic Control 23(5), 901-907.
100.2
"""No
9
10
Fig. 3. Summed squared error sequence, quality, and yield evolution in noisy tracking of EtOH, OUR, CER, and RQ . No control corresponds to 100%. against was operation with no control. From the SIMO tracking in figure 3 it is seen how the control algorithm reduces the Summed Squared Error Sequence ek (SSES) with approximately 80% in the first batch run . Furthermore, after having been trained on the first three batch runs the control algorithm reduces SSES with more than 94% in subsequent batch runs except for batch run no. 6. By coincident the direction of the perturbation in batch run no. 6 corresponds to starved yeast and is unique to batch run no. 6 and particularly intricate. However, the performance of the control algorithm in batch run no. 7 does not suffer from the performance drop (still, 40% of SSES was rejected) in batch run no. 6.
er
6. CONCLUSION
In the present paper, the framework presented by Lee et al. (2000) has been extended with a dynamic description of the batch-to-batch behavior of initial conditions and with their effect on
388
Copyright © IFAC Dynamics and Control of Process Systems, Jejudo Island, Korea, 200 I
c:
0
[>
Pu blications www.elsevier.comllocate/ifac
DEVELOPMENT OF LEARNING CONTROL FOR REPRODUCIBLE AND HIGH QUALITY OPERATION OF BATCH PROCESSES Dennis Bonne * , 1 Sten Bay
*
J~rgensen * ,2
Department of Chemical Engineering, Technical University of Denmark, DK-2BOO Lyngby, Denmark
Abstract: In this contribution, a framework for modeling and control of batch processes has been extended with the effect of initial conditions on batch evolution and with the inclusion of measurement noise. The imbedded model state estimation algorithm has been stabilized and extended with the ability to process noisy measurements. A scheme has been developed for identification of high dimensional models and it has been shown that models with immense dimensions can be obtained from sparse, noisy data. The extended framework and identification scheme were validated through simulated trajectory tracking of fed-batch fermentations . Copyright © 2001 [FAC
Keywords : Batch process, modeling, identification, model predictive control, iterative learning control, yeast fermentation
1. INTRODUCTION
dition , implementation of control will not only increase reliability, but also enable process optimization. The modeling and control framework presented in this paper is designed to reject the two first mentioned types of variation.
With an increasing industrial demand for flexible and specialized production methods, batch and semi-batch processes are becoming ever more popular. Control of batch operation, nonetheless, has received little attention in academia. Most batch processes are operated according to off-line specified operation plans, which typically specify reference trajectories for physical and chemical properties. Batch and semi-batch processes however, often deviate from normal/ desired operation due to three types of variation: Variation that occurs during a batch run, variation between batch runs, and variation that reappears after a certain period. These variations from normal operation cause operation scheduling and plant design to be conservative. Consequently, to reduce the conservatism of present operation scheduling and plant design , feedback control is required to reduce variations from normal operation. In ad1 Corresponding author; fax: +45 4593 2906 e-mail: db@kt .dtu .dk 2 fax: +45 4593 2906 and e-mail : [email protected] .dk
2. MODEL DEVELOPMENT The following batch modeling approach is based on a framework (Lee et al., 2000) where nonlinear and time-varying dynamic behavior is approximated with a set of simple, Linear Time-Invariant (LTI) models. Lee et al. (2000) suggested using a set of Finite Impulse Response (FIR) models and Gregersen (2001) suggested using a set of AutoRegresive models with eXogenous (ARX) inputs. Here the model may be set up using either FIR models or ARX models . Operation modes of batch processes will vary with time. Therefore, the number of measurements (which require trajectory tracking) and actuators may vary with time . Often it will be necessary to include additional measurements for monitoring and estimation purposes, and to measure the quality of the
and
383
the resulting combined model is:
end-product. Consequently, it seems reasonable to define the following variables:
ek+1 = ek -
• Control variable, u(t) E Rnu(t). • Controlled output variable, y(t) E Rny(t), which requires tracking control. • Secondary output variable, s(t) E Rns(t), used only for monitoring. • Quality variable, q ERn., only measured after a batch.
ek = ek
y:=[yT(l)
]T , YT (2) ... yT(N)
T
d:= [d (l)
o
ek+1
(1)
-
GY ~Uk+1 -
e~:=y-y=e~+v~,
One major drawback of the proposed modeling approach is the immense dimensionality of the resulting set of models. E.g., assuming causality a process with one output and one input sampled 100 times during a batch will be parameterized by 5050 parameters using FIR models and using ARX models 198 model orders and 594 parameters assuming all model orders set equal to 3. In practice, this immense dimensionality will render any identification problem ill-conditioned. Fortunately, there is some structural information in the resulting set of models G. Firstly, by causality the upper triangular is 0, secondly, the model parameters time variation can reasonably be assumed smooth (L 1 ), i.e., the diagonal elements in G will vary in a smooth manner, and thirdly, the impulse responses can reasonably be assumed smooth (L 2 ), i.e. , the column elements in G will vary smoothly. Regularization utilizes such a priori knowledge, reducing the amount of data needed, by introducing constraints or penalties in order to attract excessive degrees of freedom toward reasonable values. Unfortunately, regularization will result in biased estimates, i.e., there will be a tradeoff between variance and bias. The optimal tradeoff is found by selecting the regularization weights Ai that minimize an average prediction error over the validation batches. The regularized system:
+ Mni(Yini).
(3)
Grni~Yini,k+l
+ w~, (4)
Uk+1 - Uk and ~Yini,k+1 = The matrices GY and Grni are linear system approximations and e~ is the part of e~ that reappears in e~+l' Both w~ and v~ are assumed zero-mean i.i.d. (independent and identically distributed) sequences with respect to k.
where ~Uk+l
Yini,k+1 -
eZ + ~~,
3. IDENTIFICATION
The nonlinear model (3) may be approximated by a set of linear models. In the present paper N(u , d) will be approximated with one linear model for each sample instant. Let y be the specified and preferably optimal output reference trajectory and let k be the batch index. Then the output error trajectory, e Y , is defined as:
= e~
(7)
(2)
:= y(O) .
Note, not all initial conditions are measurable and/or physically meaningful. The output sequence is in general related to the input sequence, initial condition, and disturbance sequence by a nonlinear model:
e~+l
(6)
where vg is assumed a batch wise white disturbance term , z~ is the measurement of e~ , and ~~ is assumed batch wise white measurement noise.
where N is the batch length. Let the initial condition be defined as:
y = M(u, Yini , d) = N(u, d)
= ek0 + Vk0
z~ =
]T , dT (2) . .. dT(N)
Yini
+ Wk,
To enable usage of a Kalman filter for estimation of the initial condition eO, the initial condition is modelled with a simple batch wise random walk model:
The input, output and disturbance sequences over the entire batch are defined as: u T (l) ... u T (N-1) ]T ,
- G ini eZ+ 1
2.3 Model for Initial Condition
2.1 Model for Controlled Output
u:= [uT(O)
G~Uk+1
+ Vk·
=
Yini,k·
2.2 Combined Model
b reg = AregCf',
By applying the same modeling approach to the secondary and quality outputs and defining:
(8)
is solved in a least squares (LS) sense. The vector Cf' is the unknown model parameters, while Areg contains the system inputs and initial conditions (and outputs for ARX models) and breg contains
(5)
384
[!:::] = [~ ~] [!:::=~] - [~:=~] ~Uk(t Zk(t) =
y(/)
[0
Ht] [ek,t]
ek,t
- 1),
+ ~(t),
(11)
where t = 1, ... , N, Zk(t) is the measurement of the output error ek(t) at time t, and ~(t) is assumed zero-mean i.i.d. measurement noise at time t. Note, this state-space formulation does not seem observable, i.e., e.k,t is not observed in Zk(t) . However, remembering that ek = e.k+Vk, e.k ,t can in fact be observed through ek,t with a filter. At the start of a batch the states in (11) satisfy:
1(0)
Fig. 1. Block diagram of an Iterative Learning batch Control (ILC) algorithm in a Model Predictive Control (MPC) framework with combined output trajectory tracking and end-product quality control. the system outputs. The regularized system (8) is given as:
[ o~] [>.11] =
c{>,
(9)
>'2L2
Note that the structure of Ht and ek(t) depend on which measurements are available at sample t . Let €k,t+mlt be the optimal estimate of the error trajectory of the kth batch given measurements up to time t and m future control moves. Then the estimate €k ,t+mlt is given by:
where the exact system set up from Nk batch data sets:
(10)
€k,t+mlt = €k ,tlt - G'(' ~uk.t·
In literature, several approaches to selection of scalar optimal regularization weights have been presented. However, seemingly no work has been done with vectorial regularization weights based on model properties. For further discussion of how to solve the regularized LS problem (8) the reader is referred to Hansen (1996). Furthermore, usage of high quality data for the identification will decrease estimation variance. The quality of batch identification data can be determined by statistical comparison to a normal reference model obtained through multi-way analysis. For further discussion of high quality or normal batch identification data the reader is referred to Gregersen and J0rgensen (1999) and Nomikos (1996).
(13)
The optimal input movements are determined by employing the standard quadratic objective function known from conventional MPC control.
where Qt and R t are weighting matrices. The analytical solution to (14) is given by:
~uk.t = [G'('T QtG'(' + Rtt1G'('T Qt€k,tlt , (15) where only ~Uk(t) from ~uk\ is implemented. In practice input signals and their rate of change must be constrained to ensure non-provoking or system friendly and physically meaningful inputs. If possible the inputs should be constrained such that the process avoids regions of instability increasing the level of safety in the event that the controller should set out. Furthermore, the predicted outputs may need to be constrained to meet additional safety requirements and ensure operation in a physically meaningful operation region. Therefore, (14) is extended to:
4. MODEL-BASED ITERATIVE LEARNING
CONTROL Based on the newly obtained batch transition model, Iterative Learning batch Control (ILC) in a Model Predictive Control (MP C) framework can be applied for combined trajectory tracking and end-product quality control (Lee et al., 2000). Model predictive control has been comprehensively covered in literature, e.g. by Mayne et al. (2000). The block diagram of the control algorithm with the extensions introduced in the present paper (i.e., initial conditions and measurement noise) is shown in figure 1. From eq. (6) the following periodical, Linear Time-Varying (LTV) state-space model is obtained for control design:
Win
~Uk , t,Ek , t
[-2 {~uk.r Qt~uk.t + 2R; ~uk.t 1
+ELStEk,d], S.t.
At~uk.t::; €k ,t ~
385
0,
Bt(Ek,t),
(16)
'R;
where Qt = Gr;'T QtGr;' + R t and = -e.LltQtGr;' . By proper configuration of A and B both the outputs and the inputs may be constrained. The constraint softening slack variable E ensures feasible solutions and St is the corresponding weighting matrix. Since the secondary output is not controlled, the weighting matrix Qt has the following shape: Qt
=
Qi 0 [ o
1
0 0 0 0 0
Ql
Gt -
Pk,t1t-l] = [ ~k,tlt-l P k,tlt-1 P k,tlt-1
]
(20)
.m
T T
+ GiniPeoklkGini + R w ,
= Pk-1 ,NIN +GiniPeoklkG'fni + Rw
+ Rv .
In the above P is the prediction error covariance of the error trajectory, P is the prediction error covariance of the part of error trajectory that reappears, P is the cross-covariance between e and e, 'P is the covariance of the innovations, K is the Kalman gain, lC is the normalized Kalman gain, and Rn,t is the measurement noise covariance at time t . Note that the noise sequences are assumed both mutually independent and independent of the error sequences at time t . For further discussion of numerically stable implementations of Kalman filters the reader is referred to Grewal and Andrews (1993) and Kailath et al. (2000) . 4.1.1. Estimation of Initial Conditions The estimated initial condition eZ,k is obtained from: AO
ekl k =
(1
0 - "'k ) eAOk- 1Ik-1 + "'kZk>
+ Rvo) x (PeOk-1Ik-1 + Rv o + Rn o) -1, (1 - "'k) (PeOk-1Ik-1 + Rvo) ,
"'k = (PeOk-1Ik-1
PeO k1k =
(21)
where PeO k1k is the prediction error covariance of the initial condition, "'k is the Kalman gain, Rvo is the covariance of v O, and Rno is the measurement noise covariance.
1
P k,t-1 It-l]
' 5. SIMULATION RESULTS
~k,tlt] = [~k,tlt-1 ]
To investigate the applicability of the proposed modeling approach and controller design in biochemical industry a yeast producing fed-batch fermentation process was studied. To simulate fed-batch fermentation of yeast a biochemically structured model presented by Lei et al. (2001) was used. The model describes aerobic growth of Saccharomyces ceravisiae on glucose and ethanol
ek,tlt-l
-K k,dHte.k ,tlt-l - zdt)], 1t-1] lC k,t -- [Pk,t p tHT ' k,tlt-1 'Pk ,t = H t P k ,tlt-1 H ; + Rn ,t,
[gini] e2 ,k'
-
Pk,OIO
and the updating part is,
[ ek,tlt
+ Rn,t] -1 ,
= Pk-1 ,NIN + GiniPeoklkGini + Rw,
A
(18)
P k,t-1It-1 P k,t-1It-1
T
x [H t P k ,tlt-1 H ;
Pk,OIO = Pk-1 ,NIN
t:.Uk(t - 1),
[~k,t-1It-l
-
+ ICk,tKk,t'
K k,t = [P k,t 1t-1 ] H; P k,tlt-1
-
The optimal (i.e., minimum variance) estimates of the model states are obtained from a Kalman filter. Here the Kalman filter can be thought of as being made up of a set of Kalman filters or as an Extended Kalman Filter (EKF) without subsampling. By defining Rw and Rv as the covariance matrices of Wk and Vk, respectively, the predictive part of the stabilized Kalman filter (Thornton and Bierman, 1978) is,
1
-T
+Kk,tICk ,t
Pk,OIO
4.1 Kalman Filter
_ [G t -
Pk ,t1t-1] P k,tlt-1 P k,tlt-1
~k'OIO] = [~k-1'NIN] _ [ ek,OIO ek-l,NIN
The presented control technique transfers the refined control error sequence as initial condition for the next batch. Furthermore, it computes the input sequence as a change from the previous batch based on the control error prediction. This use of information from previous batches provides the control algorithm with an offset elimination action along the batch index. Thus, through this action the control algorithm is expected to eliminate offset errors both from disturbances that reappears in subsequent batches, and from model bias, i.e., the disturbances represented by w. A convergence prove for a batch MPC controller based on statespace model (11) can be found in previous work by Lee et al. (2000). This proof also holds for this extended algorithm.
[!:: :::=~] = [!:::=:::=~]
[~k,tlt-l
with the initial conditions,
(17)
.
~k,tlt Pk,t1t] = [ Pk,tlt Pk ,tlt
(19)
lCk ,t = 0.5K k,t 'Pk,t - lCk,t,
386
• The model validation was based on three normal batch data sets.
and focuses on the pyruvate and acetaldehyde branch points, where overflow metabolism occurs when the growth changes from oxidative to oxidoreductive. The fermentation was fed with substrate and the feed flow rate was manipulated as the control variable. In the simulation the ethanol concentration (EtOH), Oxygen Uptake Rate (OUR), and Carbon dioxide Evolution Rate (CER) were measured on-line, and it was assumed that there was no time delay of the measurements. From the on-line measurements OUR and CER the Respiratory Quotient (RQ) was calculated: RQ = CERjOUR. To characterize this fedbatch fermentation the following dimensionless measures were defined:
To facilitate the identification task a highly exciting input sequence was chosen. The input sequence was constructed by adding a Pseudo Random Binary Sequence (PRBS) to the nominal input trajectory. At each sample the PRBS was uniformly scaled by maximum ±1O% of the nominal trajectory. Note, highly exiting input sequences are not system friendly and thus seldom applicable in practice. For further discussion of input sequences see Pearson (1998). Since the basic assumption of this inter batch modeling approach is that trajectories vary only slightly from batch to batch, the validation batch data sets were generated with input trajectories as those used for identification, i.e., data sets similar to the identification data, but not used in the identification. From table 1 it can be seen that quite reasonable models can be obtained from noise free data for all the outputs with relative prediction fits of approximately 2 x 10- 2 . Furthermore, as indicated in table 2 models of similar quality can be obtained from noisy data (the noise level is based on laboratory experiments and industrial data) without pre-filtering. One of the three model cross-validations from the noisy identification is shown in figure 2.
• Quality == final yeast concentration minus final ethanol concentration x 1000 and substrate concentration x 200 normalized by unit concentration (1 9/l) . • Yield == produced yeast divided by fed substrate on a mass basis. The weights on ethanol and substrate concentration in the quality measure could be based on costs of downstream processing relative to the sales price of yeast, but here they are set such that variations in yeast, ethanol, and substrate concentrations are equally represented in the variation of quality. 5.1 Identification
5.2 Control
For simplicity the fermentation was modeled with EtOH, OUR, CER, and RQ as output variables and no quality variables. RQ was treated as a measurement to test if the dynamic behavior could be modelled with the proposed mode ling approach. To study the applicability of the proposed modeling approach to processes where initial conditions are not obtainable, the dynamics in G ini were discarded. A set of FIR models was identified in the following scenario:
The control simulations were designed to accentuate the iterative learning batch control feature and study the relation between trajectory control and reproducibility. The performance of a Single Input Multiple Output (SIMO) ILC algorithm without quality control was validated by randomly perturbing the initial conditions with max ±1O% in each batch run and introducing a persistent bias in the input trajectory by lowering the feed substrate concentration with 10%. The benchmark which the performance was tested
• The measurement noise was approximated by white noise with a standard deviation equal to one third of the maximum noise level. Maximum measurement noise: EtOH, ±3%; OUR,±2%;CER,±2%;RQ,±2%;u,±1%. • The initial conditions were uniformly perturbed by maximum ±1O%. • The identification was based on ten normal batch data sets. • The identification was optimized by testing possible combinations of the regularization weights ),1 and ),2 and choosing the combination that minimized the mean prediction error when cross-validated. The optimal regularization weights are shown in tables 1-2 along with their respective mean prediction fits.
Table 1. Optimal regularization weights and relative prediction fit for noise free identification. Output variable ~EtOH ~OUR
~CER ~RQ
'\1
'\2
fit
1.18 0.05 0.05 0.02
0.04 0.01 0.01 0.01
1.5746x 10 =-2" 1.9223x1O ·2 1.7818x1O =-2" 1.6298 x 10-"""7"
Table 2. Optimal regularization weights and relative prediction fit for noisy identification. Output variable ~EtOH ~OUR
~CER ~RQ
387
'\1
'\2
fit
1.18 0.05 0.05 0.26
0.05 0.01 0.01 0.02
1.5708 x 10 -2 2.4224x1O= 2.0479x 10 -2 2.3902x1O =2
batch evolution and with the inclusion of measurement noise. The modeling framework results in models with immense dimensions and hence an identification scheme has been developed utilizing structural information through regularization. Applying this identification scheme, high dimensional models can be obtained from sparse, noisy data. The model state estimation algorithm included in original framework has been stabilized and extended with the ability to process noisy measurements. The extended modeling and control framework and identification scheme were validated through simulated trajectory tracking of fed-batch fermentations with satisfactory results. Hence, the present extended framework for modeling and control of batch processes exhibits significant potential for industrial implementation.
~~] i]:::~ !j;:;~ ..,
~
o
I
2
3
..
5
•
7
e
,
10
o
I
2
3
..
5
11
7
IS
•
10
0
1
2
3
..
5
,
7
•
•
10
OO:f ~·050~---7-,--!-,-'----!3--.C-----!-'--!-,---!,:---,C---'-,---.J,O Time (h)
Fig. 2. Validation of noisy identification in batch incremented variables. True batch difference, solid line; predicted batch difference, dotted line.
i~b: : Z\ ::J !~t::52:1 !::~ ,
,
3
..
5
6
7
8
10
1
2
.,
..
5
•
7
a
10
1
2
3
•
5
6
7
a
7. REFERENCES
Gregersen, Lars (2001) . Monitoring and Fault Diagnosis of Fermentation Processes. PhD thesis . Technical University of Denmark. Gregersen, Lars and Sten Bay J0rgensen (1999). Supervision of fed-batch fermentations . Chemical Engineering 10umal75, 69- 76. Grewal, Mohinder S. and Angus P. Andrews (1993) . Kalman Filtering - Theory and Practice. Prentice Hall. N. J . Englewood Cliffs. Hansen, Per Christian (1996). Rank-Deficient and Discrete Ill-Posed Problems. Polyteknisk Forlag. Lyngby. Kailath, Thomas, Ali H. Sayed and Babak Hassibi (2000). Linear Estimation. Prentice Hall. N. J . Englewood Cliffs. Lee, Jay H., Kwang S. Lee and Won C. Kim (2000) . Model-based iterative learning control with a quadratic criterion for time-varying linear systems. Automatica 36(5), 641-657. Lei, F ., M. Rotb011 and S. B. J0rgensen (2001) . A Biochemically Structured Model for Saccharomyces cerrevisiae. Accepted for Journal of Biotechnology. Mayne, D. Q., J. B. Rawlings, C. V. Rao and P . O. M. Scokaert (2000). Constrained model predictive control: Stability and optimality. Automatica 36, 789--814. Nomikos, Paul (1996) . Detection and diagnosis of abnormal batch operations based on multi-way principal component analysis. ISA Transactions 35(3), 259-266. Pearson, Ronald K. (1998) . Input sequences for non linear modeling. Nonlinear Model Based Process Control. Proceedings of the NATO Advanced Study Institute pp. 599-62l. Thornton, C. L. and G. J . Bierman (1978) . Filtering and error analysis via the U DUT covariance factorization . IEEE Transactions on Automatic Control 23(5), 901-907.
100.2
"""No
9
10
Fig. 3. Summed squared error sequence, quality, and yield evolution in noisy tracking of EtOH, OUR, CER, and RQ . No control corresponds to 100%. against was operation with no control. From the SIMO tracking in figure 3 it is seen how the control algorithm reduces the Summed Squared Error Sequence ek (SSES) with approximately 80% in the first batch run . Furthermore, after having been trained on the first three batch runs the control algorithm reduces SSES with more than 94% in subsequent batch runs except for batch run no. 6. By coincident the direction of the perturbation in batch run no. 6 corresponds to starved yeast and is unique to batch run no. 6 and particularly intricate. However, the performance of the control algorithm in batch run no. 7 does not suffer from the performance drop (still, 40% of SSES was rejected) in batch run no. 6.
er
6. CONCLUSION
In the present paper, the framework presented by Lee et al. (2000) has been extended with a dynamic description of the batch-to-batch behavior of initial conditions and with their effect on
388