878
I n d . Eng. C h e m . Res. 1994,33, 878-888
Multivariable Iterative Extended Kalman Filter Based Adaptive Control: Case Study of Solid Substrate Fermentation John G. Sargantanis and M. Nazmul Karim' Department of Agricultural and Chemical Engineering, Colorado State University, Fort Collins, Colorado 80523
In solid substrate fermentation (SSF),the on-line measurements of the states of the fermentation like biomass content, dry matter content, and moisture content are not possible. Also, the control of the temperature and the moisture content is critical for optimization of the process. A multivariable adaptive control structure along with state estimation using an iterative extended Kalman filter (IEKF) is proposed for the control. The IEKF uses the measurements of total wet weight and C02 evolution rate to estimate the states. An autoregressive with exogenous inputs (ARX) model is used as a model of the process relating controlled variables with manipulated variables (dry air flow rate and the water replenishment rate) and forms the basis for the determination of future control inputs. The simulation results show that a better control of the moisture content can be achieved when compared to the single input-single output (SISO) control strategy. Introduction
In solid substrate fermentation (SSF) the microbial growth and product formation occur on the surface of the solid material. The limitations of SSF are, among others, the inhomogeneous conditions in the reactor and the difficulty of heat removal due to the low thermal conductivity of the fermentedmatter. SSF has the advantages of superior productivity, low waste water output, and improved product recovery (Hesseltine, 1977, 1987). Several conventional methods have been used to remove excessheat (Naraharaet al., 1984;Silman, 1980). Of them the most promising is the method of evaporative cooling (Ryoo, 1990) using partially saturated air. The cooling is effected by water evaporation from the substrate to the air stream. Various models have been used for the growth and product formation in SSF. Also, modeling of the heat effects has been reported (Saucedo-Castaheda et al., 1990). In this study, the logistic model of growth (Mitchell et al., 1991) without any product formation (SSF for protein formation) is employed. For this particular system, which will be described later on, work has been done on the kinetics (Barstow, 1987), the modeling (Ryoo, 1990; Sargantanis et al., 1993),and the control of the temperature (Ryoo, 1990)by an adaptive controller. A modification of this scheme is proposed in this work to simulataneously control the temperature and the moisture content. The control strategy consists of two steps. The first is the open-loop model building based on input-output mapping of the system using ARX (autoregressive with exogenous inputs (manipulated variables)). Other parametrized models that have been used include NARMAX models (Billings and Fadzil, 1985; Chen et al., 1990b; Haber and Unbehauen, 1990) (N stands for nonlinear and ARMAX for autoregressive moving average with exogenous outputs), artificialneural networks (ANN) (Chen et al., 1990a), and radial basis function networks (RBF). NARMAX and RBF can deal with severe nonlinearities of the process. The model parameters are identified using a recursive least squares algorithm (Ljung, 1987). Issues in recursive off- or on-line identification and on model selection have been extensively covered in the literature (Diekmann and Unbehauen, 1979; Gupta and Fairmann, 1972;Ljung and Soderstrom, 1983; Sagara and Zhara, 1988). The second
* To whom all correspondence should be addressed. E-mail:
[email protected]. OSSS-50S5/94/2633-0878$04.5~I0
step is to close the loop, design and implement an adaptive controller based on the minimization of an objective function containing a weighted s u m of squares of the errors of the controlled variables over a prediction horizon and a penality term for the future manipulated inputs. The MPC (model predictive control) is one of the many different adaptive control schemes that have been used (i.e., model reference control, self-tuning regulator, selfoscillating control). Also, states in the SSF are estimated using an iterative extended Kalman filter, which is used in conjunction with multivariable ARX identification and adaptive control. The iterative Kalman filter is expected to reduce the error associated with the approximation of a nonlinear process by first-order derivatives. The simultaneous state estimation with adaptive control has not yet been attempted for the SSF process. State estimation has been combined for other processes with model predictive control (Ricker, 1990). Also, a stateestimation based MPC has been developed by Lee et al. (Lee and Morari, 1992). Process Description
The microorganism used for the SSF was Rhizopus oligosporus (NRRL 2710). Flour-free yellow corn grit (no. 25 mesh particle size) was used as the substrate. The purpose is to produce cell protein, hence we want to maximize cell growth. Approximate composition of the substrate wasprotein,7.0%;fat,4.9%;starch,78.2%;fiber, 1.3%;ash, 1.7% w/w on a dry basis. The mineral medium used had the followingcomposition (in g/100 mL): (NH4)2SO4,7.59; urea, 4.09; NaHzP04,1.726; MgSOg7H20,0.059; KC1,0.159; CaC12,0.05%FeSOg7H20,0.137;alanine, 0.045. The preparation of the material for the fermentation and the inoculation procedure are described in (Ryoo, 1990). The fermentation vessel was a rocking reactor made of acrylic plastic. The reactor had the following dimensions: inner diameter, 5.1 cm; outer diameter, 12.7 cm; length, 12.2 cm; and holding volume, 1.296 L. The rocking vessel has inside and outside voids to provide uniform water supply and aeration. Air and water were supplied from the inner void part through fine holes. Gravity causes the water to permeate through the substrate. The reactor was rotating in a clockwise-counterclockwise sequence, with 1 revolution per 5 min. The temperature of the incubator and the inlet air temperature were controlled by independent on-off loops, using an 0 1994 American Chemical Society
Ind. Eng. Chem. Res., Vol. 33, No. 4,1994 879
Figure 1. Experimental setup: 1, reactor; 2, motor; 3, weighing plate; 4,balance; 5, incubator; 6, water bath; 7, humidity sensor; 8, heat exchangers; 9, water condenser; 10, drying column; 11, COz analyzer; 12, drying column, 13-15, humidifying columns; 16, air filter; 17-18, mass flow controllers; 19, personal computer; 20, thermocouples.
electrical heater along with a fan and a heated water bath, respectively. Total air flow rate and air relative humidity were controlled by blending two streams, one of which was totally humidified and the other dried by passing through a silica gel column. The two mass flow meters measuring the dry and wet air flow rates were regulated by a computer. The total air flow was controlled as the sum of the dry and the wet air flow rate, and the relative humidity was regulated by the dry air flow rate. The water replenishment flow is regulated by a small pump, and its rate was to be determined by a controller based on an ARX (autoregressive with exogeneous inputs) model of the process. The computer used an analog card for reading various temperatures in the reactor, incubator, air inlet, and outlet, air moisture at the outlet, and flow rates of dry and wet air. The total wet weight of the fermenter was measured by a balance, and the C02 evolution rate was monitored through an infrared COZanalyzer. They were also read on-line by the computer. The computer can also give control signals to the manipulated variables. The experimental setup is shown in Figure 1.
Model The dynamic behavior of the SSF process is described by unsteady-state differential equations for the following variables: biomass content ( X I ) ,total dry matter including biomass (xz), moisture content (XS),temperature of the fermented matter ( x 4 , total biomass weight (xg),total wet weight @e), and C02 evolution rate (CER). For the cell growth a logistic type of model was used (Okazaki and Sugama, 1979; Okazaki et al., 1980; Mitchell et al., 1991). The cell growth is described by two parameters, the maximum specific growth rate ~ L Mand the maximum , no substrate inhibition attainable biomass content x ~with term included. For the dry matter (substrate and biomass), the substrate consumption for growth and maintenance was taken into account. The parameters Y S / Xand m are determined using simple stoichiometric formulas for the anabolism and catabolism, approximating the substrate as consisting of starch and NH3 as the nitrogen source. The molecular formula of the dry biomass was CH1.781Oo.636N0.079 (Alvarez-Martinez, 1987). The equation for the variation of moisture content includes the water production for maintenance and growth
as well as the water exchange with the inlet air. Steadystate air flow and saturation of the air at the exit temperature are assumed. The activity coefficient of the water of the fermented matter is taken to be 1. The exiting air is also assumed to be in thermal equilibrium with the fermented matter at any time. The thermal balance takes into account the evaporative heat loss, the heat transfer with the air in the incubator, and the heat production associated with growth and maintenance. The equations for the total biomass weight and the total wet weight are easily formulated from the equations for the biomass content, total dry matter, and moisture content by using the relationships x5 = ~ 1 x 2and x6 = xz(1 + x 3 ) . The equation for the CO2 evolution rate has terms for the growth and maintenance production of C02. In all equations, homogeneity in the reactor (substrate composition, water content, biomass content, and temperature) is assumed. The model equations are:
dx = (hA(x, - !Pew) + F(E, - E,) dt
+
The yield and maintenance coefficienta for the substrate consumption, and water, COz, and heat production are taken from Ryoo (1990). The total air flow rate fl was kept content at 11L/min. The total pressure was kept at 676.2 mmHg, and the inlet air temperature was kept at 37 "C. The calculation of the absolute humidity of the air was done by multiplying the relative humidity of the air with the saturation humidity calculated using an Antoine equation. For the inlet air, the relative humidity rhin is given by fl-u, rhin = fl The enthalpy of the air contains contributions for the sensible heat of the wet air and the latent heat of
880 Ind. Eng. Chem. Res., Vol. 33, No. 4, 1994 Table 1. Values of the Constant Model Parameters parameter value units g of dry matter/g of biomass -0,59989 g of dry matter/g of biomaswh -0.01932 g of water/g of biomass 0.304 g of water/g of biomaswh 0.0106 kcal/g of biomass 1.9098 0.0787 kcal/g of biomaswh g of COdg of biomass 0.76394 g of COdg of biomaswh 0.031 ;;"i L/min 11.0 hA kcallh.' C 8.451 OC 37.0 37.0 OC C S 0.00036 kcal/g-OC C S 0.0010 kcal/g'C
evaporation of the water. The variation of C(M and with temperature and moisture content is given by:
XM
where subscripts T and W denote temperature and moisture dependence, respectively. The expressions for CCMW,PMT, XMW, and XMT have been developed (Sargantanis et al., 1993) and are given by pMW
CCMT=
exp(0.447 - 5 . 8 1 0 ~- ~5.284~: - 1.556~:) exp(-2.242
Then the estimated values of the controlled variables at time step k are
where 8 is the estimated parameter vector and y ( k ) is the vector of the controlled variables. The parameters of the linear model 8 can be estimated by various methods. We use the recursive least squares identification with exponential forgetting factor (Astram and Wittenmark, 1989; Ljung and Soderstrdm, 1983)to account for the dynamics of the process (batch operation). Before we set on the recursive identification of the parameter vector, an initial estimation of the parameter vector is obtained according to the followingprocedure. Initially the followingmatrices are formed N j=l N
(11)
+ 0 . 0 5 2 ~- ~0.772 x 10-3~,24.002
xMW
water flow rate, respectively). Of the two controlled variables, temperature is readily measured, while moisture will be estimated using an iterative extended Kalman filter. Also let us define the parameter vector B(k):
X
10-7x,3) (12)
= exp(4.143 - 20.776X3 + 1 8 . 9 6 0 ~ 4.906~3~) ~~ (13)
= exp(266.986 - 2 2 . 5 7 5 ~+~0.630~2- 0.006~,3) (14) For moisture contents greater than 66% (in wet basis), thevaluesofpMWandxMwaretakentobe constant. Finally, the heat capacity of the fermented matter is calculated as a weighted average of the heat capacity of dry matter and the heat capacity of the water: xMT
cp =
+ CP*X3 1 x3
CP,
+
The values of the various constant parameters are given in Table 1.
where X is the exponential forgetting factor, which serves to give less weight to the old data points. The initial time stepj = 1corresponds to the time to + max(n,m)At, where to is the initial time and At is the sampling time. This gives enough initial sampling intervals to form the first 4 vector, 4(1) = @(to + max(n,m)At). Then after N time intervals, the initial estimates of and the identification error covariance matrix P(N) are calculated from
P(N) = v,(N)-' &N) = V p 9 - l
v,y(N)
(21)
(22)
After the intial estimate is obtained, the parameter matrix and the-covariance matrix are updated for the time step k given B(k-11, P(k-1), and the output vector y ( k ) by the following equations:
e(k) = y(k) - P(k-1) 4 ( k ) u2 =
Adaptive Control
X2
+ &k)
P(k-1) $(k)
(23) (24)
The process described above can be modeled using a linear autoregressive with exogenousinputs (ARX) model. In this model the values of the controlled outputs at a time step k are related to the n previous values of the outputs and the previous m values of the manipulated variables. Let us define the 2(n + m) vectors:
d(k) = [Yl,k-l,
.**)
Yl,k-nt Y2,k-1,
Y2Jt-nt U1&l, UlJt-m, UZ,k-l,
*.*?
...) U2,k-ml
T
(16)
where y 1 and y 2 are the controlled variables (temperature and moisture content, respectively) and u1 and u 2 are the manipulated variables (dry air flow rate and replenishment
After identifying the parameter matrix, the manipulated variables can be determined for the current time step by appropriately minimizing an objective function. Let J be NP
Nc-1
J=1
t=O
J = z e ( k+ j ) TI',e(k+j) +
AuT(k+i) I',Au(lz+i)
(27)
Ind. Eng. Chem. Res., Vol. 33, No. 4, 1994 881 I' and R are 2 N p X 2Np and 2 N p X 1 matrices given where e(k+j) = y(k+j) - r(k+j), u(k+i) = [UlC+i, ~ ~ h + i l ~ where , Au(k+i) = u(k+i) - u(k+i-1), and r(k+j) is the set point by: vector at the time step k +j. N p is the prediction horizon, Jrlo and Nc is the controller horizon. Constant control u(k+Ncl) = u(k+Np-1) = ... = u(k+Np-1) is implemented for the (Np- Nc) future control inputs after time tk+N,-1. Np was chosen to be 3 while Nc was selected to be 1. This means that only the control input u(k) that is going to be implemented for the next time interval needs The solution to the problem that uses the inputs u instead to be calculated. The objective function is a weighted of their change Au is given by the same expression except sum of the deviation of the controlled variables from the for the omission of the term rzu(k-1). For Np = 1, rl= set point and contains a contribution from the future Iand I'2 = 0, the controller is reduced to an one step ahead control inputs. The matrices l'l and r2 are diagonal and controller. The procedure that was followed to run the are used as a weighting factor between the accuracy of adaptive control simulation was (i) For an initial time of control and the magnitude of the control input. They can 1.5 h, an initialization was performed and the matrices be tuned to yield the desired performance of the controller. V4(N)and V,(N) were built. At time t = 1.5 h the initial The outputs at future time steps are estimated using the error covariance matrix P and 4 were formed. (ii) Until parameter matrix available at the current time step (certainty equivalence approach). This contains a sigtime t = 6 h open-loop identification was run. This is nificant amount of uncertainty, so that the extrapolation necessary for the estimates of the parameter matrix to to distant time intervals is likely to be inaccurate. For NC converge. For the first two steps, random independent equal to 1 and constant future control up to time tk+ivp1, perturbations of the manipulated inputs were applied. using sequentially the relation y(k+j) = &%)[y(k+j-l), The inputs had mean values of 3 L/min and 6 g/h, y(k+j-2), u(k+j-1), u(k+j-2)]T for j = 1, ..., Np, the respectively, and had a probability of 50% to vary in the expression for the second-order ARX model (m= n = 2) next time interval. The level of variation was for both is: *20% of the maximum value. (iii) After 6 h, the loop was closed and the manipulated inputs were determined by Y(k) = the control law and implemented for the next time interval. The sampling time was 0.1 h, and the value of the forgetting factor X was 0.995. = A&) + A2y(k-l) + A,u(k-1) + Bu(k)
-
-01
. I
Iterative Extended Kalman Filter (28)
where Y (k)is an extended output vector over the prediction horizon with dimensions 2Np X 1. The matrices AI, At, As, and B (allwith dimenions 2 N p X 2) can be expressed as functions of the partitioned matrix &k) = [818283841T (Np = 3), where the matrices Bi have dimensions 2 X 2:
e,(e,e, + e,) + 419,+ 0, +
+ 8,
The purpose of the Kalman filter is to estimate states of the system which are difficult to measure experimentally, using readily available measurementa. Also,reduction of the noise corruption of the measurements can be achieved using a Kalman filter. In SSF,cell content, dry weight, and moisture content cannot be measured on-line. Available experimental measurements for the particular system are (1) the carbon dioxide evolution rate, 21, and (2) the total wet weight, 22, through a balance. Estimates for the biomass content x1 and the total dry matter x 2 can be obtained by the above two measurements 21 and 22. Moisture content can be indirectly estimated from the combination of the filtered measurement 2 2 and the estimated dry matter weight x2. The discrete time extended Kalman filter is a minimum variance estimator which minimizes the variance between the states estimated by a continuous nonlinear.mode1and the states estimated from the discrete-time nonlinear equations of the observer. The model can be written in the following form (Lewis, 1986)
1
The objective function to minimize J can be written compactly as J = (Y- R)TT(Y - R) + AuT(k) I'2Au(k). Then the solution to this unconstrained minimization problem is
~ ( k =) [BTI'B+ r2]-'(BTI'(R - A,y(k) - A,y(k-1) Asu(k-l) + I'2~(k-l))) (33)
where x(t) = [ x l ( t ) , x 2 ( t ) l T is the 2 X 1 state vector, f(x(t),u(t)) = [ f l , f i l T is the 2 X 1 state function vector, z ( t 3 = [21,22ITisthe 2 X 1 observation vector at discrete time t k , g(x(tk),u(tk-l)) = [g1,g2IT is the 2 X 1 observer equation vector, w ( t )is the state model error vector (2 X l),and ~ ( t kis) the measurement noise vector (2 X 1). The elements of the f(x(t), u(t)) and Z ( t k ) vectors are
882 Ind. Eng. Chem. Res., Vol. 33, No. 4,1994
(37)
where weh is the initial wet weight of the fermented matter. The extended Kalman filter algorithm follows the following steps (i) Initialization: P(0) = Po,x(0) = ~0 (expected value of x at time t = 0) (ii) Prediction:
where
df5,.. =
-(ax,a -1 -)ddtx ~
=
xs
x(klk-1)
x(k-l(k-1)
+ Sa_lf(x(k-llk-l),u(k-l))
dt (41)
+ Skk_JA(x(k-llk-l)) X + P(k-1Jk-1) AT(x(k-llk-l)) + Q(k))dt (42)
(1- (Ys/x)(x,))2
P(klk-1) = P(k-llk-1) P(k-llk-1)
For the state prediction we integrate using a second-order Runge-Kutta. (iii) Correction:
G(k)' = P(k(k-1) HT(x'(klk)) [H(x'(k(k)) P(k(k-1) X HT(xf(klk))+ R(k)l-' (43) x'+'(k(k) = x(k(k-1) + G ( k ) ' ( ~ (-k g(x'(k(k),u(k-l)) ) + H(x'(klk))(x'(klk) - x(klk-1))) (44)
P(klk) = [I - G(k)' H(x(k(k))lP(klk-1)
(45)
where P(kJk)is the a posteriori estimation error covariance matrix 2 x 2 at discrete time step k, Q and R are the 2 X 2 covariance matrices of the model and measurement errors, respectively. G ( k ) is the Kalman gain matrix at time step k, A(x( )) is the Jacobian matrix of the state equations and H(x( )) is the Jacobiaamatrixof the observer equations (both 2 X 2). The elements of the Jacobians of the states and the measurements are
(47)
The correction includes an iterative procedure (Bellgardt et al., 1986) which s t a r b with xl(k(k) = x(klk-1) and proceeds until (x'+l(k(k)- x'(lz(k)(Ic, where c is a vector of small numbers. At convergence, the covariance matrix P(k-llk-1) is corrected according to eq 26 using the converged Kalman gain matrix G(k). The noniterative filter is just a special case of the above algorithm, in which eqs 24 and 25 are applied just once and they are not allowed to converge. After the measurement vector is filtered, an indirect estimate of the moisture content is obtained through the expression gl(x(k(k1,U(tk-1)) = (1+ x&lk))x2(klk), which can be solved for xdklk)
x,(klk)
gl(x(klk), u(k-1)) -1 x,(klk)
(54)
Results and Discussion The two-input two-output adaptive identification and control and the estimation of the states in the SSF are implemented at the same time. The flow diagram of the implemented scheme is shown in Figure 2. The simulation is done off-line. The model described earlier, eqs 1-7, is used to simulate the process. The integration is done with a Runge-Kutta fourth-order method with a step size of 0.001h. The outputs of the process, the temperature, and the two measurements are corrupted with random white noise. The level of the noise is f0.4for both measurements. The two measurements, along with the two manipulated inputs delayed by a time step, enter the iterative extended Kalman filter routine. Estimates of the cell content, the
Ind. Eng. Chem. Res., Vol. 33, No. 4, 1994 883
I
2nd order ARX
I
5 1
I
I
1.1 1 2 ..... 2:1 . . '
I 4 - i
PROCESS
2.2 3.1 - - 3.2 -.. 4.1 .. .
4.2 _.....
3 - $1
I eI I + '
ITERATIVE EXTENDED KALMAN FILTER
Figure 3. Parameters of the identificationfor a second-order ARX model, first four rows of 8. 2nd order ARX
2,
I
1.5
Figure 2. Proposed structure for the adaptive control of temperature and moisture content in SSF combined with Kalman filter estimation ofthestates of theprocess;ltl,&,.&are theestimated biomaascontent, dry matter, and moisture content, respectively,and x g and x7 are the simulated C02 evolution rate and total wet weight, respectively.
I
0.5 0
dry matter, and the moisture content are obtained. The current temperature and estimated moisture content enter the ARX identification step with inputs and outputs delayed by one and two time steps, respectively. The identification produces estimations of the parameter matrix, which is fed to the controller, along with the desired set point values of the controlled variables. The set points can be set arbitrarily, as will be done in this study, or they can be determined by an optimization routine that can use the Kalman filter estimates of the cell mass in an attempt to maximize the cell production over a period of time. The determined manipulated inputs are fed back to the model and the Kalman filter to complete the loop. All parts in the loop are interrelated. The time step employed is the same for both identification and the Kalman filter. The diagonal elements of the matrices in the Kalman filter are initialized to P(1,l) = 0.0005,P(2,2) = 0.0001, Q(1,l) = 0.002, Q(2,2) = 0.0001, R(1,l) = 0.40, and R(2,2) = 0.40. The values of the parameters PM and X M in the Kalman filter simulation are assumed to be constant, because the Kalman filter assumes constant mean values ( p =~0.21 and X M = 0.49) in order to test the ability of the Kalman filter to give correct estimates of the states, when some parameters are wrongly calculated. The manipulated inputs cannot be negative and have maximum values of 11 L/min (total air flow rate) and 20 g/h (replenishment water flow rate), respectively. The optimization for the calculation of the manipulated variables is unconstrained, but when the calculated values exceed the limits, the manipulated variables are reset to the corresponding limit, with a small white noise superimposed to drive the identification. The implemented set point changes are as follows: the temperature is set to 37 "C for the first 25 h, then to 38 "C from 25 to 38 h, and then back to 37 "C, while the moisture content set point is 60% for the first 20 h, then 61 % from 20 to 33 h, and finally 60%. The simulation was run for 50 h. The variables in both identification and the Kalman filter were normalized around some mean values to ensure comparable magnitudes of the variables in the multivariable routines. The basic results of the simulation are as follows:
-0.5
-2
'
I
IO
5
0
15
20
25 Time (h)
30
35
40
4.5
SO
Figure 4. Parameters of the identificationfor a second-order ARX model, last four rows of 8. MIMOZndARX,IEKF,Np;.3,h l = l , h24.5,h3.0. h 4 4 , lam4.995 0.3
-u "x
't
0.2
0.1
0
I -0.1 -0.2
-0.3
0
5
10
IS
20
25
Time (h)
30
35
40
45
50
Figure 5. Identification error of temperature with a second-order ARX model, Np = 3, iterative Kalman filter, rl = diag(l.0,0.5) and r2 = 0.
The identification performs very well. Although it is unstable at the beginning of the open loop with the random values of the manipulated variables, the parameters converge and behave smoothly in the closed loop. The values of the parameters B are given in Figures 3 and 4 for 50 h and second-order ARX. The identification time is short, as the process is a t unsteady state (batch operation) and it needs to be controlled at an early stage. The identification error of both temperature and moisture content is shown in Figures 5 and 6, respectively. Except at the moment when the loop is closed both errors show good behavior. The identification error of the temperature
884 Ind. Eng. Chem. Res., Vol. 33, No. 4,1994
p 3 E
E .5
MIMO 2nd ARX. IEKF, Np=3. h l = l . hZ.0.5. h3.0, h4.0. lam.0.995
MIMO 2nd ARX. IEKF, Np=3. h l = l . hZ.0.5. h3.0. h4.0. lamdl995
0.05 0.04
I
1401
I
0.03
0.02
120 100
0.01
80
I
,
V
I
1
,
I
I
i
I
,
kalman process -----
-
0
60-
4.01
4.02
40-
4.m
1
20 0
4.05
0
5
10
I5
20
25 30 Time (h)
35
40
45
SO
Figure 6. Identification error of moisture content with a secondorder ARX model. Np = 3, iterative Kalman filter, rl= diag(l.0.0.5) and l2' = 0.
0
0.6 .___._. ------
5
IO
15
20
25 30 T i e (h)
35
40
45
SO
Figure 7. Estimation of biomass content X I by iterative extended Kalman filter.
15
20
25 30 Time (h)
35
40
45
SO
MIMO 2nd ARX. INT, N p l . hl-1. hZ-O.5, h3-0, h4-0. h - O . 9 9 5
38.5
1
1
36l 0
0
IO
Figure 9. Estimation of the total biomass weight x 3 by iterative extended Kalman filter. 39
MIMO 2ndARX.IEKF. Np=3, h l = l . h2.0.5. h 3 4 , h4=0,lam=0.99.5
5
I
scl poult
'
5
'
10
' IS
' 20
'
'
25 30 Time (h)
'
'
35
40
' 45
50
Figure 10. Temperature ("C). Control withone step ahead controller with rl=diag(l.0,0.5),I'z = 0, and Kalman filter forced to very small gain (integration). MIMO 2nd A m . m . N p - 1 . hl-I. h2.0.5. h3=0. h M . la1n-0~995
MIMO 2nd ARX, IEKF. N-3. h l = l , hZ.0.5. h3.0. h4=0.lam.0.995
0
-
m
0
5
10
I5
20
25 30 Time (h)
35
40
..._ --..__ 45
50
Figure 8. Estimation of the dry matter x2 by iterative extended Kalman filter.
is larger, but the future prediction of the temperature is better with the second-order ARX model than the future prediction of the moisture content. The identification was hardly affected by the variation of the forgetting factor A. The iterative extended Kalman filter did not perform as well as expected. It is noisy at 10 h of fermentation for the estimation of the dry matter (Figure 8). The biomass content is predicted better, but the estimation is biased (Figure 7). The estimation of the dry biomass weight is good (Figure 9). It seems like small uncertainties in the biomass content that are caused by the inaccurate PM and XM amplify the error in the dry matter estimation, which in turn makes the indirect estimation of the moisture content behave poorly also. Possible reasons for this behavior of the Kalman filter can be (1)the great variation of the growth parameters PM and XM with the temperature
5
IO
15
20
25 30 Time (h)
35
40
45
50
Figure 11. Dry air flow rate, u1, (L/min). Control with one step ahead controller with I'l = diag(l.0,0.5), z'I = 0, and Kalman fiiter forced to very small gain (integration).
and moisture content, (2) the errors introduced by the prediction step in the iterative extended Kalman filter. The integration is much less accurate in the Kalman filter than in the process model integration, and (3) there is interaction between the Kalman filter and the identification. The iterative Kalman filter (converged in about 10 iterations) was better than the noniterative Kalman filter. The influence of the initial choice of the error covariance matrix P,the model covariance matrix Q,and the measurement covariance matrix R is too complicated to be analyzed. Choosing small P and Q we force smaller Kalman gain and we reduce the correction and the noise in estimation. The same uncertainty holds with thechoice of the mean growth parameters. The estimation is very sensitive to p~ and XM. The EKF is known (Lewis, 1986) to diverge when an appreciable modeling error exists. The IEKF reduces the divergence. Possible actions against modeling uncertainties and colored noise are measurement
Ind. Eng. Chem. Res., Vol. 33, No. 4, 1994 886 MlMOZnd~INT.Npl,hl-l.W.S.h3d).h4-0.LM1d.995 0.68
PP
0.64 0.62
1
........
,
MIMOZnd AFX, INT.Np-I. hl-1.0. hZ-0.2. h3-0.02, h4-0. lun-0.995
12,
1
1
,
I
,
,
,
,
*
,
............
. ....................... ...._.....
-
....
056
05*
I 3 ; 1
054 0 0.52
pmss
set point
.--.-
ka--
5
10
15
24l
25
30
35
40
45
-2
5(1
0
5
10
I5
20
Time (h)
Figure 12. Moisture content (% wet basis). MIMO control based on Kalman filter estimate. Control with one step ahead controller with rl = diag(l.0,0.6), rz = 0, and Kalman filter forced to verysmall gain (integration).
25
30
35
40
Time (h)
45
50
Figure 15. Dry air flow rate, ul,(L/min). Control with one step ahead controller with I'l = diag(l.0,0.2), = diag(0.02,0.0), and Kalman filter forced to very small gain (integration). MIMOZndARX.INT.Npll,h1-1.0,h2d.2,h3-0.02,h4-0.lund).995 0.68
0.66
'S f
0.54
'
-5L
0
5
'
IO
' 15 .
' 20
' 25
Time (h)
' 30
' 35
' 40
'
45
Figure 13. Replenishment water flow rate, 112, (glh). Control with one step ahead controller with rl=diag(l.0,0.5), rz = 0, and Kalman fiter forced to very small gain (integration). -0
2nd A N .
INT.N o d . hl-1.0.
' 5
' 10
'
I5
'
20
' 25
'
30
' 35
'
40
' 45
IO
5
'
"
I5
20
"
25
Time (h)
30
35
"
40
45
50
Figure 16. Moisture content (5% wet basis). MIMO control bawd on Kalman filter estimation. Control with one step ahead controller with r1= diag(l.0,0.2), rz = diag(0.02,0.0), and Kalman filter forced to very small gain (integration). MIMO2ndARX. INT.Np-I. hl-1.0, h2-O.2. h3-O.02. h4-0. iam-O.995
W . 2 . h3-0.02. h4-0. iam-0.995
Time (h)
"
0
25
36I 0
1
O.X1
'
50
..........
0.62
'
511
5
5 1 0
"
5
IO
"
I5
20
' 25
"
30
35
"
40
45
I 50
Time (h)
Figure 14. Temperature (OC). Controlwithoneetepaheadcontroller with rl=diag(l.0,0.2), = diag(0.02,0.0),and Kalman filter forced to very small gain (integration).
Figure 17. Replenishment water flow rate, uz, (g/h). Control with one step ahead controller with rl=diag(l.0,0.2),2'l = diag(0.02,0.0), and Kalman filter forced to very small gain (integration).
noise excitation of the states, exponential data weighting, and expansion of the state vector using a colored noise filter. The error of the model is likely to be colored, but it is difficult to simulate a colored noise model in such a dynamic and changing environment because the error of the model depends on the control moves through the different trajectories of temperature and moisturecontent. Thus, forcing the Kalman filter to small values improves the short-time behavior (local up and downs) of the filter but cannot reduce the mismatch that is caused by the uncertainty of the model parameters. The control performance is strongly dependent on the K h a n filter estimation of the moisture content. A singleinput single-output temperature controller (with the second manipulated variable at a constant mean value with a superimposednoise) performs very well. The noise in the temperature measurement is taken to be f0.2 "C. Though this noise, f0.2 "C,seems high, in the real-time
experiments,the nonhomogeneity of the fermented matter causes these types of fluctuations in the temperature measurement. The temperature and moisture content of the process are highly coupled. An increase in the moisture content increasesthe heat capacity of the fermentedmatter and thus deteriorates the heat removal capacity. Also, if the dry air flow rate increases, the evaporation causes moisture removal, which drives the second manipulated variable to increase too. The control of temperature in a SISO manner is feasible. The purpose of the MIMO approach is to try to control the moisture content, too. Achieving good control of the temperature in the presence of moisture content control and maintenance of the moisture content within prescribed limits would be a satisfactory result. By employing a SISO adaptive control approach with no Kalman filter and on-off control of the moisture based on a mass balance of water in the fermentor with an estimation of the rate of water production from
886 Ind. Eng. Chem. Res., Vol. 33, No. 4, 1994
,-.
u
3 738 5
MIMO 2nd ARX. IEKF, Np.3, hl=l, hZ4.5. h 3 4 , h 4 4 l a m 4 . 9 9 5
MIMO2nd ARX.NCT, INT.Np-3. hl-I, h2-0.2. h 3 4 . h46).lani=0.995
38.5
12
1
I
1
i
i
7
36.5
36
1 1 , 0
5
,
,
,
,
10
15
20
25
, 30
TIME (h)
, 35
, 40
,
I
45
50
Figure 18. Temperature ("(2). Control with three step ahead (NP = 3) controller with rl = diag(l.0,0.2), r2 = 0, and Kalman filter forced to very small gain (integration). 10,
MIMO2nd ARX.NCT, IHT, Np=3, h l = l . h2=0.2, h 3 4 , h 4 4 . l a m 4 . 9 9 5 , , I , I
I
I
I
I
,
0
5
10
I5
20
ZS
Time (h)
30
35
.
IO
IS
20
25
30
35
40
45
setpoint
MIMO 2nd ARX, IEKF, Np=3, h l = l , h 2 4 5 , h 3 4 , h4=0,larn=O 995
36
o
s
IO
15
zo
zs
Time (h)
30
3s
40
4s
,
-
process .---
SO
TIME (h)
Figure 19. Dry air flow rate, u1, (L/min). Control with three step ahead (Np = 3) controller with r l =diag(l.0,0.2), rz = 0, and Kalman filter forced to very small gain (integration).
SO
MIMO 2nd ARX,IEKF,Np=3, h l = l . h24.5. h 3 4 . h4=0,lam4,99S 6 8 1 , , , , , , , ,
khan
5
45
Figure 21. Dry air flow rate, u1, (L/min). Control with three step ahead (Np= 3) controller with rl = diag(l.0,0.5), r2 = 0, and iterative extended K h a n filter. 0
0
40
0
5
10
I5
20
25
30
Time (h)
35
'
40
4.5
SO
Figure 22. Moisture content (7% wet basis). Estimation baaed on Kalman filter estimation. Control with three step ahead (Np = 3) Controller with rl = diag(l.0,0.5), I'2 = 0, and iterative extended Kalman filter.
so
Figure 20. Temperature ("C). Control with three step ahead (NP = 3) controller with F l = diag(l.0,0.5), r2 = 0, and iterative extended Kalman filter.
the COz evolution rate yielded an increasing with time moisture content profile with a final value of around 70% in wet weight basis (Figure 24). Higher moisture content results in lower cell biomass production, as the optimum moisture content lies between 60 % and 66 % (Ryoo, 1990). The multivariable approach succeeds in better containing the variation of moisture content. The control of the temperature is far better than the control of the moisture content. Both are impaired by the somewhat biased estimates of the Kalman filter. The effect of the prediction horizon, Np, in the control suggests that larger Np improves the control of temperature by rejecting moisture content set point changes (Figure 18) but impaires the moisture content control, due to the poor prediction of the ARX model caused by the noisy estimation of the moisture content by the Kalman filter. The coupling of the inputs is apparently not well estimated
0
5
10
IS
u)
25 30 Tune (h)
35
40
45
SO
Figure 23. Replenishment water flow rate, u2 (g/h). Control with three step ahead (Np = 3) controller with rl= diag(l.0,0.5), r2 = 0, and iterative extended Kalman filter.
for long prediction horizons. An increased order ARX identification or the employment of NARX models (Chen et al., 1989; Chen et al., 1990b) would possibly alleviate these problems. The manipulated variables u1 and u2 are inversely influenced by increasing N p . The first shows less variation and becomes more sluggish, while the variation of the second is increased. By changing the elements of I?l we can put more weight on either temperature or moisture content control. By properly tuning the matrix rz the controller can reject the disturbances in temperature control affected by a change in the moisture content set point (compare Figures 10 and 14). The elements of I'g need to assume small values especially the one associated withug, which is usually set to 0. Generally,
,
Ind. Eng. Chem. Res., Vol. 33, No. 4, 1994 887 Nomenclature
SlSO controlled SSF. i d k t estimation of mois. coni.
biomass content (g of biomass/g of dry matter) total dry matter (g) x3: moisture content (g of water/g of dry matter) x4: temperature of the fermented matter ("C) x5: biomass weight (9) (= ~ 1 x 2 ) total wet weight (9) (= xz(1 + x 3 ) ) CER: COZevolution rate (g/h) u1: flow rate of dry air portion (L/min) UZ: replenishment water flow rate (g/h) p ~ maximum : biomass specific growth rate (h-1) XM: maximum attainable biomass content (g of biomass/g of dry matter) I&,: absolutehumidity of the inlet air (g of water/g of dry air) H,,,,absolute : humidity of the outlet air (g of water/g of dry air) &: enthalpy of the inlet air (kcal/g of dry air) Eout:enthalpy of the outlet air (kcal/g of dry air) Y S I Xyield : coefficient (g of dry matter/g of biomass) m: maintenance coefficient (g of dry matter/g of biomassqh) Y H ~ oyield : coefficient for water (g of water/g of biomass) mn,?: maintenance coefficient for water (g of water/g of biomass-h) YQIX:energy yield coefficient (kcal/g of biomass) mQ: energy maintenance coefficient (kcal/g of biomass-h) YCO; yield coefficient for COZ(g of COdg of biomass) rnco,: maintenance coefficient for C02 (g of COz/g of biomass-h) F: total dry air flow rate (g of air/h) fl total air flow rate (L/min) rhin: relative humidity of inlet air hA: overall heat-transfer coefficient (kcal/h-"C) Taw:temperature in the incubator ("C) Tinlet: inlet air temperature ("C) cp: average heat capacity of the fermented matter (kcal/g"C) cp,: heat capacity of the dry matter (kcal/g"C) cp,: heat capacity of the water (kcal/g°C) web: initial wet weight of the fermented matter (g) XI:
x2:
3
0.56 0.58 0.54
i1
0.SZL
o
"
s
io
'
IS
"
20
zs
Time (h)
1
'
30
35
"
40
4s
1 I
so
Figure 24. Moisture content ('7% wet weight) for SSF with SISO adaptive temperature control and indirect estimation and on-off control of the moisture content with the set point at 60% moisture content.
the control of temperature is made better by increasing Np or by tuning r2,while the moisture content is better controlled by a step-ahead controller. The choice depends on the objectives set. If temperature control is important a step-ahead controller with an adjustment in the matrix r2 is the most proper. Finally, a larger noise in the Kalman filter predictions makes interaction phenomena larger especiallybetter the moisture content and the temperature control. The control using the values of the inputs u in the objective function is better than the one that uses the change of the inputs Au, and this is the one implemented in all the results shown. Figures 10-13 show the temperature, dry air flow rate, moisture content, and water replenishment flow rate as determined by a one step (Np = 1) ahead controller (r2 = 0) with diagonal elements hl and h2 for rl. Figures 14-17 show the same variables determined by a controller with Np = 1 and diagonal elements hl and h2 for rl and h3 and hq for r2. Figures 18 and 19 show the temperature and dry air flow rate calculated with the previous matrix settings but with Np = 3. All these controllers have been calculated with a simple integration in the Kalman filter in order to show the performance of control alone. Figures 20-23 show the temperature, dry air flow rate, moisture content, and water replenishment flow rate with control in conjunction with the Kalman filter. The parameters of the controller are the same as those in Figures 10-13. Finally, Figure 24 shows the moisture content variation by using a SISO control of the temperature and indirect estimation and on-off control of the moisture content with the set point at 60% moisture content. Conclusions The following final conclusions can be drawn: An ARX model can be implemented in SSF to serve as a basis of a MIMO control of temperature and moisture content. The control of both variables is shown to be feasible. The K h a n filter can estimate well the biomass content, but the uncertainty of the growth parameters forces a noisy estimation of the moisture content. The containment of the moisture content within some limits is better with the MIMO control than by on-off controlling the moisture independently of the temperature control. The process is shown to be coupled, and the effects of both temperature and moisture content need to be taken into account. Improvement of the Kalman filter prediction would likely lead to better control and optimization of the process.
Literature Cited Alvarez-Martinez, L. Modeling fungal growth on extruded corn by Solid Substrate Fermentation, PhD thesis, Colorado State University, 1987. Astrom, K. J.; Wittenmark, B. Adaptive control; Addison-Wesley: Reading, MA, 1989. Barstow, L. Evaporative temperature and moisture control in Solid Substrate Fermentation. Master's thesis, Colorado State University, 1987. Bellgardt, K. H.; Kuhlmann, W.; Meyer, H. D.; Schuegerl,K.; Thoma, M. Application of an extended Kalman filter for state estimation of a yeast fermentation. ZEE Proc. 1986,133(51,226-234. Billings, S.A.; Fadzil, M. B. The practical identification of systems with nonlinearities. In Proceedings of the 7th ZFAC Symposium on Identification and System Parameter Estimation; Pergamon Press: Oxford, England, 1985;pp 155-160. Chen, S.; Billings, S. A.; Luo, W. Orthogonal least-squares methods and their application to non-linear system identification. Znt. Control 1989,50 (5),1873-1896. Chen, S.;Billings,Grant, P. M. Nonlinear system identification using neural networks. Znt. J. Control 1990a,51,1191-1214. Chen, S.;Billings, S. A.; Cowan, C. F. N.; Grant, P. M. Practical identification of NARMAX models using radial basis functions. Znt. J . Control 1990b,52 (6),1327-1350. Diekmann, K.;Unbehauen, M. Recursive identification of multiinput multi-output systems, In Proceedings of the ZFAC Symposium, Darmstadt; Pergamon Press: Oxford, 1979;pp 423-429. Gupta, R. D.; Fairman, F. W. Periodic sampling and continuoustime system identification. Znt. J. Sci. 1972,3,187-196. Haber, R.; Unbehauen, H. Structure identification of nonlinear dynamic systems-A survey on inputloutput approaches. Automatika 1990,26,651-677. Hesseltine, C. W. Solid Substrate Fermentation: An Overview.Znt. Biodeterior. 1977,23,79-89.
888 Ind. Eng. Chem. Res., Vol. 33, No. 4, 1994 Hesseltine, C. W. Microbiology of oriental fermented foods. Ann. Rev. Microbiol. 1987, 37, 575-601. Lee, J. H.; Gelormino, M. S.; Morari, M. Model Predictive Control of Multi-Rate Sampled-Data Systems: A State-Space Approach. Znt. J. Control 1992,55, 153-191. Lewis, F. L. Optimal estimation: with an introduction to stochastic control theory; Wiley: New York, 1986. Ljung, L. System Identification, Theory for the user; Prentice H a k Englewood Cliffs, NJ, 1987. Ljung, L.; SBderstrBm,T. Theory and practice of recursive identification; The MIT Press: Cambridge, MA, 1983. Mitchell, D. A.; Greenfield, P. F.; Doelle, H. W. An Empirical Model of Growth of Rhizopus oligosporus in Solid-state Fermentation. J. Ferment. Bioeng. 1991, 72 (3), 224-226. Narahara, H.; Koyama,Y.; Ywhida, T.;Atthaeampunna, P.; Taguchi, H. Control of water content in a Solid State culture of Aspergillus oryzae. J. Ferment. Technol. 1984,62,463-459. Okazaki, N.; Sugama, 5.Growth estimation of Aspergillus oryzae cultured in solid media. J.Ferment. Technol. 1979,57,408-412. Okazaki,N.; Sugama, S.;Tanaka,T. Mathematical model for Surface Culture of Koji Mold. J . Ferment. Technol. 1980,58,471-476. Ricker, N. L. Model predictive control with state estimation. Znd. Eng. Chem. Res. 1990,29 (3), 374-382.
Ryoo, D. On-line estimation and control in Solid Substrate Fermentation. PhD thesis, Colorado State University, 1990. Sagara, S.; Zhara, Z. Y.Numerical integration approach to on-line identification of continuous system in presence of measurement noise, In Preprints of the 8th ZFAC Symposium on Identification and System Parameter Estimation; Pergamon Press: Oxford, England, 1988; pp 27-31. Sargantanis, J. G.; Karim, M. N.; Murphy, V. G.; Ryoo, D. Effect of operating conditionson Solid Substrate Fermentation. Biotechnol. Bioeng. 1993,42 (2), 149-158. Saucedo-Caatafleda,G.; Gutierrez-Rojas,M.; Bacquet, G.; Raimbault, M.; Vinieqa-Gonzalez, G. Heat Transfer Simulation in Solid Substrate Fermentation. Biotechnol. Bioeng. 1990,35 (E), 802808.
Silman, R. W. Enzyme formation during Solid Substrate Fermentation in rotating vessels. Biotechnol. Bioeng. 1980,22,411-420. Receiued for review April 15, 1993 Revised manuscript received October 22,1993 Accepted December 21,1993. a Abstract
15, 1994.
published in Aduance ACS Abstracts, February
