Ind. Eng. Chem. Res. 1992,31, 848-855
848
Tedder, D. W.; Rudd, D. F. Parametric studies in industrial distillation, Part 11. Heuristic optimization. AIChE J. 1978b,24,316. Thompson, R. W.; King, C. J. Systematic synthesis of separation schemes. AIChE J . 1972,18,941. Umeda, T.; Niida, K.; Shiroko, K. A thermodynamic approach to heat integration in distillation systems. AIChE J. 1979,25,423. Underwood, A. J. V. Fractional distillation of multicomponent mixtures-Calculation of minimum reflux ratio. J. Inst. Pet. 1946,32,614.
Underwood, A. J. V. Fractional distillation of multicomponent mixtures. Chem. Eng. h o g . 1948,44(8),603. Westerberg, A. W.; Stephanopoulos, G. Studies in process synthesis-I. Branch and bound strategy with list techniques for the synthesis of separation schemes. Chem. Eng. Sci. 1975,30, 963.
Received for review November 4, 1991 Accepted November 14,1991
Dynamic Modeling and Optimal State Estimation Using Extended Kalman Filter for a Kraft Pulping Digested Chimmiri Venkateswarlu and Kota Gangiah* Simulation, Optimization and Control Group, Indian Institute of Chemical Technology, Hyderabad 500 007, India
A nonlinear dynamic model of delignification has been developed for kraft pulping of Pinus casia wood, based on 36 batch experiments at various operating conditions. Further, the validity of the mathematical model is ascertained with independent experiments. Extended Kalman filter and nonlinear observer are applied to estimate lignin and carbohydrate contents during pulping by employing the developed model in conjunction with the effective alkali concentration as available process measurement. The dynamic simulation results indicate that extended Kalman fii$er is more efficient than the nonlinear observer for estimating known and unknown s t a b which enable obtaining the desired pulp quality. Introduction Batch digesters are used to cook wood chips with a chemical solution called liquor to chemically dissolve the organic lignin binder to free the individual cellulose fibers. In the kraft process, the lignin bond is made soluble by cooking the chips in an alkali-based liquor, which release the wood fibers. The alkali liquor contains two active chemicals, sodium hydroxide and sodium sulfide. The degree of delignification of a pulp is generally measured by an indirect laboratory analytical test known as the Kappa number. Direct determination of lignin and carbohydrates through analysis during pulping is tedious and time-consuming. Thus,it is difficult to control the product quality by using only infrequent and delayed laboratory data, whereas liquor samples during pulping can be drawn easily and can be analyzed for alkali concentration. In kraft pulping, the desired pulp quality is achieved by monitoring the degree of delignification. Batch digesters are subject to disturbances such as cooking liquor changes, production schedule changes, poor heat distribution, changes in the initial starting temperature, and changes in the length of cooking time. In the past, many studies of the kraft pulping process have been made to develop theoretical and experimental models of varying complexity for control and design purposes. A state space model which considers time-varying pulping process variables is quite useful for process design and operation. Further, this kind of model can cope up with the stochastic features of the process such as modelling errors, sampling uncertainties, and process noise. The operating behavior of the pulping digester can be predicted by using the available process measurements in conjunction with the mathematical model of the process. Certain techniques known as estimators can be employed to filter the process measurements and to estimate the
* Author to whom correspondence should be addressed. t IICT Communication
No. 2731.
lignin content, carbohydrate content, and alkali concentration with time during a cook. The variables thus estimated can be used for timely monitoring the pulping process. In this study, the objective is to develop a dynamic model for the kraft pulping digester and to apply the extended Kalman filter and the nonlinear observer to estimate lignin and carbohydrate contents by using effective alkali concentration as available process measurement. The performances of the state estimators are evaluated through simulation studies. Experimental Aspects Experiments are carried out to generate data which can be useful to construct a rigorous dynamic model of a batch digester. The laboratory digester of M/K Systems, Inc., U.S.A., which is shown in Figure 1, consists of a pressure vessel, a liquor circulating pump, a heat exchanger, a pressure gauge, and a temperature indicator. Pinus casia logs from Manipur forests of India are made into smaller size chips of about 25-mm length, and the screened sample of 3-5-mm width is used as a raw material. The wood is analyzed (Tappi Standards and Suggested Methods, 1962) for lignin, carbohydrate, and alcohol-benzene solubility as 27.90%, 67.90% and 4.37%, respectively, on an oven dry basis. A woad sample of 250 g (moisture-free)is taken for each experiment and is added to the circulating liquor of specified concentration with a wood to liquor ratio of 1:4. Thirty-sixbatch experiments are carried out at three temperature levels (160,170, and 180 "C),four effective alkali concentration levels (35,40,45,and 50 g/L) at each temperature and three cooking time levels in the range of 15-270 min at each alkali concentration. In each experiment, cooking temperature is specified in the temperature range of 160-180 "C for a specified period of time. The specified temperature is controlled in terms of a set point by using on-off control. Hence, temperature rises from an initial value, say 30 "C, to the set-point temperature
0888-5885/92/2631-0848$03.00/00 1992 American Chemical Society
Ind. Eng. Chem. Res., Vol. 31, No. 3, 1992 849
x
Figure 1. Schematic of laboratory digester: 1, digester with cover; 2, basket; 3, insulation;4, thermistor probe; 5, electric heater; 6, flow control valve; 7, pressure gauge; 8, temperature indicator; 9, pump; 10, variable transformer; 11, top.
and pulping continues at that point until the specified period of time. White liquor is prepared by using 110 g/L as 96% NaOH and 65 g/L as 50% Na2S in 3:l ratio. The liquor is analyzed by using a pH meter titration procedure (Macdonald, 1969). The liquor composition is Na2S as Na20 = 27.59 g/L, NaOH as Na20 = 82.305 g/L, effective alkali as NazO = 96.1 g/L, and sulfidity = 25.1%. White liquor of specified concentration on an oven-dry wood basis is charged to the digester in each experiment. A wood sample is added to the circulating liquor of specified concentration, and the continuous heater is immediately switched on. Liquor circulation is carried continuously through a circulating pump. The temperature is observed, and the effective alkali concentration is determined by drawing cooking liquor samples. Liquor samples which are drawn during pulping are analyzed for effective alkali concentration by using a pH meter titration procedure (Macdonald, 1969). At the end of the experiment, pulp is removed from the digester, washed with water, and analyzed for lignin and carbohydrates (Tappi Standards and Suggested Methods, 1962). The data of carbohydrates and lignin of the wood sample at the beginning of the experiment and carbohydrates and lignin of the pulp sample at the end of the experiment are available. For each experiment, cooking liquor composition and temperature for every 5 min during the rising period, and for every 15 min during the remaining cooking period, are also obtained.
Dynamic Modeling The development of computers and statistical analysis of experimental data has made it possible to easily fit a regression equation to collected cooking data (Hatton and Keays, 1970; Hinricks, 1967). However, linearized versions of higher order regression models did not find much use as these models are not based on kinetics of the pulping proms. Experimental models are those which assume that the pulping reaction rates are kinetically controlled. The H-factor model (Vroom, 1957) is one of the earliest experimental models in which the extent of delignification in kraft pulping is expressed with time and temperature. The H-factor model is more suitable, if all the charging
conditions are constant and delignification is assumed to follow zero-order kinetics. Normally charging conditions are not constant and can vary substantially from cook to cook. In order to satisfy these conditions, other experimental models (Hatton, 1973; Chari, 1973; Venkateswarlu et al., 1986) have attempted to correlate initial charge conditions with the H-factor for predicting end-point Conditions. These nonrigorous models can be quite useful over a narrow range of operating conditions and are suitable for specific applications. Various theoretical and experimental models have been reviewed by Jutila (1980). A theoretical mathematical model consisting of a series of partial differential equations describing the combined diffusion and kinetics within a wood chip during kraft pulping has been derived and solved using the orthogonal collocation method (Gustafson et al., 1983). Model predictions are compared with experimental data reported in the literature. Such a generalized, detailed theoretical model may not be very useful in practice. However, if diffusion effects within a wood chip are not accounted for and only kinetics are considered, then the distributed parameter model proposed by Gustafson et al. (1983) reduces to lumped parameter model. In this study, therefore, a lumped parameter or state space model development for kraft pulping is attempted since this kind of model is well suited for operation and control. The experimental data of the previous section are used for the development of the dynamic model of delignification, which is similar to the state space model given by Jutila (1980). The model is described by the following equations: dC/dt = -ql exp(b - a/T)CL; C(0) = Co (1) dL/dt = -q2 exp(b - a/T)CL; L(0) = Lo (2) dCO/dt = q 3 dL/dt; CO(0) = COO (3) where C is effective alkali concentration, L is lignin concentration, CO is carbohydrate concentration, and T is absolute temperature. b and a are specific constants of a wood which resemble Arhenious equation type constants such as frequency factor and energy of activation. Optimum coefficients of the pulping model ql, q? qs, b, and a are computed for each experiment by using the Nelder and Mead simplex method (Kuester and Mize, 1973) with the following objective function:
where C is the observed value of the effective alkali concentration, C is the computed value of the effective alkali concentration, Lf is the observed end value of the lignin content, Lf is the computed end value of the lignin content, COf is the observed end value of the carbohydrate content, and is the computed end value of the carbohydrate content. The Nelder and Mead simplex method has the advantage of being able to "searchn by tentatively expanding or contracting vertices of the simplex. For the present model, the simplex has six vertices. The computed values used in the objective function are obtained by solving the pulping model equations numerically with the help of the Runge-Kutta fourth-order technique. The initial values of model Coefficients are the same for all the experimentsand are given by q1 = 1.0 x 10-5 q2 = 1.0 x 10-5 43 = 1.0
mf
850 Ind. Eng. Chem. Res., Vol. 31, No. 3, 1992 Table I. Optimum Model Coefficients expt no. 91 1 2.00189 x 10-5 1.83916 X 2 1.88171 X 3 2.00089 x 10-5 4 1.83003 X 5 1.96670 X 6 2.27134 X 7 1.65962 X 8 2.17627 X 9 2.44153 X 10 2.36954 X 11 1.79121 X 12 1.84476 X 13 14 1.66064 X 1.47889 X 15 1.40497 X 16 17 1.99552 X 2.01762 X 18 2.24152 X 19 1.52395 X 20 21 1.79417 X 22 1.57736 X 1.70573 X 23 1.66261 X 24 1.83309 X 25 1.67229 X 26 1.77417 X 27 1.62396 X low5 28 1.59936 X 29 2.01011 x 10-5 30 2.21761 X loT5 31 1.55002 X 10" 32 1.56967 X 33 1.66972 X 34 1.44311 X lo" 35 1.57606 X 36
Qz 1.70544 X 1.52003 X 1.72480 X 1.84363 X 1.58277 X 1.80996 X 1.64355 X 1.55855 X 1.58838 X 1.88117 X 1.67439 X 1.53541 X 1.51568 X 1.42037 X 1.36494 X 1.71441 X 1.81972 X 1.81229 X 2.01772 X 1.33083 X 1.73993 X 1.51945 X 1.62198 X 1.72132 X 1.50543 X 1.69669 X 1.72992 X 1.51229 X 1.56221 X 1.81992 X 1.64040 X 1.60449 X 1.34674 X 1.77002 X 1.50353 X 1.60672 x
Table 11. Statistics of Model Coefficients coeff mean std dev 91 1.82432 X 2.69876 X lo4 92 1.63235 X 1.53063 X lo4 43 1.13161 8.3801112 X b 44.5537 0.1516688 a 16 113.187 0.354 592 14
Since a and b values of Pinus casia are not known, reported values of spruce wood (Vroom, 1957) as a = 16113.0 b = 43.2 are used as initial values for the optimization routine. The optimum values of model coefficients of 36 experiments are shown in Table I. The mean and standard deviation of the model coefficients of 36 experiments are shown in Table 11. The optimum values of model coefficients of individual experiments of Table I are closely scattered with their corresponding mean values. The a and b values obtained in Table I1 for the Pinus casia wood are not much Werent from the a and b values of spruce wood (Vroom, 1957). The mean values of model Coefficients can be reasonably employed which satisfy all the batch experiments. Therefore, the dynamic pulping model which is developed in this study is described by eqs 1-3. Estimators The need for better understanding of a process is advancing the use of new estimation techniques. Although these techniques seem most often used for on-line data processing, they can also be used for off-line estimation. Filtering of measurements is also part of the estimation process, which eliminates any random noise and yields
93
10"
loT5
10" 10-5
1.237 29 1.17966 1.10354 1.24319 1.04106 0.981 82 1.269 27 1.17325 1.147 93 1.275 31 1.34562 1.21831 1.14162 1.18201 1.05473 1.14050 1.00100 1.10222 1.11221 1.10142 0.991 69 1.09902 1.024 02 1.102 29 1.10697 1.11001 1.10001 1.099 91 1.10121 1.21001 1.277 37 1.11001 1.102 13 1.11001 1.10183 1.10576 PROCESS NOISE
I"
b 44.6320 44.8164 44.7205 44.6671 44.7480 44.6712 44.595 1 44.5868 44.3946 44.6016 44.4355 44.5066 44.6423 44.4778 44.5504 44.4861 44.7818 44.4812 44.6442 44.3717 44.8525 44.6436 44.5122 44.4822 44.4501 44.4001 44.7525 44.6012 44.5442 44.4221 44.6882 44.3910 44.1459 44.3889 44.4789 44.3729
a
16 113.320 16 113.189 16 113.017 16 113.301 16 113.230 16 113.113 16 112.858 16 113.575 16 113.360 16 113.750 16 113.221 16 113.214 16 113.334 16 113.210 16 113.215 16 113.734 16 113.301 16 113.470 16 113.684 16 113.283 16 113.226 16 113.502 16 113.075 16 113.101 16113.121 16 112.476 16 111.310 16 112.796 16 113.031 16 113.480 16 113.793 16 113.578 16 113.241 16 113.238 16 113.443 16113.011
OBSERVAT ION NOISE
R'
KALMAN FILTER OR
NONLINEAR
t-
-U OBSERVER
Figure 2. State estimation scheme for the kraft pulping digester.
smooth and reliable estimates. Estimation techniques are widely used in chemical and biochemical engineering for state estimation and parameter identification (Caminal et al., 1987; Jo and Bankoff, 1976; Ramirez, 1987; Wells, 1971; Stephanopoulos and San, 1984) or process fault detection and diagnosis (Park and Himmelblau, 1983; Dalle Molle and Himmelblau, 1987). In this study, the well-known extended Kalman filter (Caminal et al., 1987) and the recently proposed nonlinear observer described in the Appendix (Frank, 1987) are applied for state estimation by using the dynamic pulping process model which is presented in the previous section. The scheme of state estimation is shown in Figure 2. The mathematical model of the pulping digester can be expressed by the following state space form X ( t ) = f ( x ( t ) , t )+ w ( t ) , x(0) = x, (5) where w is the vector of the process noise and x, is the initial condition of the state x. The linear measurement relation is given by Y(tk) = ElX(tk) + V ( t A (6) The nonlinear measurement equation can be written as Y(tk) = h ( X ( t k ) ) + V ( t k ) (7)
Ind. Eng. Chem. Res., Vol. 31, No. 3, 1992 851 where h is a nonlinear function of state x and v is the vector of observation noise. The state vector, x(t), of eq 5 can be estimated from the known process measurements, y ( t k ) ,of eq 6 with the help of nonlinear estimation techniques. The statistical expectations of w(t), v(tk)and x(0) are with the following relations: E[w(t)] = 0 E[x(O)l = xo
0
+ V(tk)
(8)
(9)
As carbohydrate is a function of lignin (eq 3),the extended Kalman filter is implemented by considering the state vector as (X)T
= (CL)T
(10)
The dynamic model of the pulping digester is used to represent f(x) of eq 5. The measured variable, y, of eq 9 is used to estimate C and L. The state estimate CO can be obtained from the estimated L as eq 3 is dependent on eq 2. A first-order nonlinear observer (Appendix) is adopted for state estimation of the pulping digester. Instead of using for local estimation, the observer is applied to estimate the states simultaneously. The available state x, and the unavailable state xu (eq A.l) can be represented by
(aT
(X,)T
=
(X")T
= (LIT
I
2.0
1
3.0
I 6.0
Figure 3. Comparison of independent experimental lignin values ( 0 )with thoee predicted by the proposed model (-) and eq 14 (---).
=
The matrices Po,Q(t),and R(tk)are used to reflect errors in the initial state, process model, and process measurements. Po, Q(t),and R(tk)are generally selected as design parameters, and Q ( t ) and R(tk)can be referred to as Q and R. Effective alkali concentration, C, with noise is used as known process measurement y: Y ( t k ) = C(tk)
I
1.0
T i m e , hours
~~(xo-X(O))(xo-X(O~)l= Po
E[v(tk)vT(tJl = R(tk)
.= 20 ._ B 10
E[V(tk)l = 0
E[W(t)WT(t)l
50
(11) (12)
Observer computations are carried out by considering the state vector given by eqs 11 and 12. The state carbohydrate content is computed by using the estimated lignin content. The matrix K, of eq A.17 is computed as
Results and Discussion A state space model involving lignin and carbohydrate content can be found to give better applicability for monitoring the process, if these variables are made available during cooking. It is possible to estimate these variables with time during a cook with the support of a dynamic model, by using suitable estimation techniques. Hence, 36 batch experiments are carried out to generate data to develop a dynamic model of delignification. The coefficients of the model, eqs 1-3, are determined by using the Nelder and Mead simplex method. In order to compare the results of model of thisstudy, another well-known model of delignification which is the basis of the H-factor
approach is also presented. This model is based on the Arrhenius-type equation, according to which the influence of pulping temperature on the rate of delignification is expressed by -dL/dt = exp(b - a / T ) (14) The optimum coefficients a and b of eq 14 are computed for each experiment by using the Nelder and Mead simplex method with a = 16113.0 and b = 43.2 as initial values which are also used in the earlier modeling. The objective function employed is given by J = (Lf - Lf)2 (15) The mean values of 36 experiments, a = 16 113.437 with standard deviation 0.3164149 and b = 39.1426 with standard deviation 0.1182407, are employed as model coefficients. The models thus developed are used for simulation to reflect the actual operating conditions. Initial effective alkali concentration and temperature are chosen as main operating variables. The following variations are introduced for the operating Variables to ascertain the effect of each variable on the digester states: effective alkali concentration = 35-50 g/L; maximum temperature = 160-180 OC. The temperature of the digester rises from the initial value to a maximum specified value (2'-) for a particular period of time Oh). In actual practice, temperature versus time history will be known. However, in the simulation environment the temperature versus time relationship is obtained by fitting to the 18 equations in two parameters (Venkateswarlu et al., 1986). The best fitted one is given by T = { (32024.58t + 1038.82)'' T, fat>t,i,
for 0 < t < tM
(16)
T is temperature in degrees Celsius and t is time in hours. The tti, is set by the specified T-. The dynamic models of eqs 1-3 and eq 14 are solved with the fourth-order Runge-Kutta method to generate true process state trajectories. The numerical simulation results of the model of this study and the other model described by eq 14 are compared with the actual independent experimental values of lignin. Since there is no provision of drawing pulp samples from the digester during a cook, independent experiments are carried out by discontinuing a cook of 160 "C,45 g/L initial effective alkali concentration, and 4 h cooking time at different cooking times starting with the same initial and operating conditions to obtain pulp samples for lignin analysis. The experimental values of lignin are closer to the results of the proposed model as shown in Figure 3. The reason for not satisfying the results of model eq 14 with the experimental lignin values may be due to employment of only accessible end lignin values of each experiment to determine the coefficients a and b.
852 Ind. Eng. Chem. Res., Vol. 31, No. 3,1992 M r
20 0
eo Y)
5z
w
I
c
40
bz
o
E o w 4 0 8 0
5 I < >
I
r Carbohydrate (*I. on wood )
1
60
E 4 0 c < @ 20 0 80
Co rbohydrate (% on wood)
(d)
60
40 20 ('1. on wood
0 0 TIME ,hours
0.5
1.0 T I M E , hours
2.0
1.5
Figure 4. Estimated states by extended Kalman filter for 160 O C range when initial effective alkali concentrations are (a) 35, (b) 40, (c) 45, and (d) 50 g/L; ( 0 )experimental, (-1 true, (- - -) estimated.
Figure 5. Estimated states by nonlinear observer for 180 O C range when initial effective alkali concentrations are (a) 35, (b) 40, (c) 45, and (d) 50 g/L; ( 0 )experimental, (-) true, (---) estimated.
The model of eq 14 is simple, but not suitable for estimation as it is not a function of a dependent variable to compute filter or observer coefficient matrices. Moreover, lignin measurement is difficult and time-consuming during pulping, whereas measurements with time are needed for estimation. The model of this study is more suitable for estimating known and unknown variables during pulping by using easily available alkali measurement. Hence, the proposed model is employed to estimate the known and unknown states for different values of operating conditions by applying both the extended K h a n filter and nonlinear observer. Random Gaussian numbers (Shannon, 1975) are used in the form of random Gaussian noise. A zero mean random Gaussian noise with a standard deviation of 0.75 is used to corrupt the true values of effective alkali to obtain alkali measurement and also to corrupt the true temperature. The values of noise covariance matrices selected for the filter are given by
The efficiency of the extended Kalman filter and nonlinear observer are analyzed by using a performance criterion which is a sum of the squared errors of true and estimated states as given by
0.5
=[ o
0
OS]
R = 0.5
(19) The design parameter of the observer K, is selected as
I =
N
c (C(t,) - C ( t k ) ) 2 + (L(tJ - i ( t k ) ) 2 + k=l ( W t k )
- C O ( W (21)
where C and 6 are true an4 estimated values of effective alkali concentration,L and L are true and estimated values of lignin content, and CO and CO are true and estimated values of carbohydrate content. The noise covariance matrices of the filter and the design parameter of the observer which are given above are selected to obtain the same performance value at the following base case: initial effective alkali concentration = 35 g/L; cooking temperature = 160 "C;cooking time = 4 h. Initial values of lignin and carbohydrate contents used in the simulation are experimentally determined values for the wood chips. Estimated states by the extended Kalman filter for various initial values of effective alkali concentration at 160 "C are shown in Figure 4. Estimated states by the nonlinear observer for various initial values of effective alkali concentration at 180 "C are shown in Figure 5. The results show that the estimated states by both the extended K h a n filter and nonlinear observer are in good agreement with their corresponding true values. In Figures 4 and 5, only the measured effective alkali is compared with the simulation or estimation results. Estimated states are filtered values of known effective alkali and unknown
80
!gz Y z
Ind. Eng. Chem. Res., Vol. 31, No. 3, 1992 853
,
I
... 40M 20
0' 0
I
1
I
1.o
I
2.0
3 .O
4 .O
;;@of
I
60
i.-
20
A
A
-
0
0
3z 80
1.0
2 .o
3 .O
4.0
1.0
2 .o
3 .o
4.0
0
0
1 .o
2 .o
0
1.0
2 .o Time ,hours
3.0
L .o
3.0
4.0
I
0
-5 U
0
T i m e , hours
Figure 6. Comparison of states for 160 OC range and 45 g/L initial effective alkali concentration: ( 0 )independent experiments, (-) true, (- -) estimated by extended Kalman filter.
-
lignin and carbohydrate with time during pulping. Further, independent experimental results (not the ones used to determine model coefficients) are used to compare the measured and unmeasured estimation variables. A random Gaussian noise with a standard deviation of 0.75 is also employed to reflect errors in the process model by corrupting temperature and model equations. A zero mean random Gaussian noise with a standard deviation of 1.0 is also used to corrupt the true values of effective alkali to obtain alkali measurement. The estimated states by the extended Kalman filter and nonlinear observer for a case of 160 O C , 45 g/L initial effective alkali, and 4-h cooking time are shown in Figures 6 and 7. The results show that the estimated lignin and carbohydrate are in reasonable agreement with their corresponding actual experimental values. At least the initial comparison shows that the extended Kalman filter is slightly better than the observer. Although the experimental values of effective alkali are on lower side initially, the latter results nearly followed the estimates. The performances of the filter and observer are evaluated interms of a quantitative performance criteria, eq 21. The performances of the estimators are judged by selecting the design parameters (eqs 17-19 for filter and eq 20 for observer) so as to provide the same performance criteria at the base case and then by evaluating the performance criteria at various operating conditions. The performances of filter and observer, when initial alkali concentration varies from 35 to 50 g/L for 160 O C (4-h cooking time), 170 "C (3-h cooking time), and 180 O C (2-h cooking time) are shown in Figure 8. The performance of extended Kalman filter is better than that of the observer for 160 "C. Although the performances of both the estimators are almost the same at the base case, the extended Kalman filter gave a slightly higher performance at 170 and 180 "C with the same initial alkali concentration as that of the base case. The trend in performance indices of the extended Kalman fiter for 170 and 180 "C is different from that of the base case, whereas the observer has shown a similar trend for all the temperatures. The process and observation noise
Figure 7. Comparison of states for 160 OC range and 45 g/L initial effective alkali concentration: ( 0 )independent experiments, (-) true, (- - -) estimated by nonlinear observer.
25
50 25 0
t
(a) I
I
II
F 35
1
1
40
45
(C) 50
Eftective 01 kali conccnlmtion(g / I )
Figure 8. Performances of extended Kalman filter (EKF) and nonlinear observer (NLO)at variow initial alkali concentrations: (a) 160, (b) 170, and (c) 180 "C.
covariance matrices, Q and R ,are the same values as in the base case. Hence, the different trend in filter performance can be analyzed in terms of the effect of operating temperature and initial effective alkali concentration. When temperature rises to the higher side from the base case, its effect on the predicted filter covariance matrix is causing a rise in error and performance index, whereas rise in initial effective alkali concentration in the presence of rise in temperature is giving weight to Kalman gain through predicted covariance and measurement Jacobian matrix which gives a slightly faster convergence and low-
864 Ind. Eng. Chem. Res., Vol. 31, No. 3, 1992
ering performance index. Since there is no involvement of noise covariance matrices in the observer computation, this type of behavior is not found in the performance indices of the observer. The differences in performance indices of the extended Kalman filter between the starting and final values of effective alkali are lower for the three temperature ranges when compared to that of the observer. The extended Kalman filter is therefore preferred for state estimation in the kraft pulping digester.
Conclusions Thirty-six experiments are carried out at various operating conditions for the kraft pulping of Pinus casia. A nonlinear dynamic model, which accounts for effective alkali, lignin, and carbohydrate concentrations has been developed by using experimental data and by optimally determining the model coefficients with the help of the Nelder and Mead optimization technique. Dynamic simulations are carried out to estimate the known and unknown states of the pulping system by employing two different methods, known as the extended Kalman filter and the nonlinear observer, and by using effective alkali as the available process measurement. The performances of the filter and observer are evaluated by selecting the design parameters to provide the same performance values at the base case and then by varying the initial effective alkali concentration. The extended Kalman filter has performed slightly better than the nonlinear observer and can be applied effectively for unknown state estimation in the pulping digester. Nomenclature C = effective alkali concentration (g/L) CO = carbohydrate content (% on oven-dry wood) Cor= end value of carbohydrate content (% on oven-dry wood) f = nonlinear function of states H = measurement Jacobian matrix h = nonlinear measurement equation Z = performance criteria used for estimators J = objective function used for the optimization routine K, = design parameter used in the observer L = lignin content (% on oven-dry wood) Lf= end value of lignin content (% on oven-dry wood) T = absolute temperature (K) t = time v = observation noise vector w = process noise vector x = state vector x, = available state vector xu= unavailable state vector y = measurement vector Appendix. The Nonlinear Observer The state vector x in the nonlinear process model can be expressed in terms of available state by measurement x, and unavailable state by measurement xuas given by ( x ) =~ (XaXJT (A.1) The nonlinear process model, eq 5, and the measurement model, eq 7, can be expressed in the following form: x = f(x,,x,) + w (A.2)
with xathe n, X 1 state vector, xuthe nu X 1 state vector, ya the ma X 1 output vector, and yu the mu X 1 output vector unavailable by measurement. f, C, and C, are considered partially differentiable nonlinear functions. The observer is described by 4 = f(&%J + Ka(ya-9,) + KU(YI;-QJ (A.7) Ya = Ca(ia) (A.8) Y u = Cu(fu) (A.9) Xa(0) = i a O (A.lO) XJO) = i u o (A.ll) and (i)T
=
(iaiu)T
(A.12)
The estimated state vector i, of fr is used to replace yu. The estimation error equation which is linearized around the state estimates 2, and 2, is given by e, = F,e, + Fueu (A.13) e, = xu - xu (A.14) A
(g 2)l
(A.16) Fu = - KU *as" The stability of the observer may be assured by an appropriate choice of the feedback matrix K,. The matrix K, is in general nonconstant and can be chosen as K u = -df -/ x-1dC,# (A.17) dxu *e,*" dxu *,,*"
where the # sign indicates pseudoinverse.
Literature Cited Caminal, G.; Lafuente, J.; Lopoz-Santin, J.; Poch, M.; Sola, C. Application of Extended Kalman Filter to Identification of Enzymatic Deactivation. Biotechnol. Bioeng. 1987,29,366-369. Chari, N. C. S. Integrated Control System Approach for Batch Digester Control. Tappi 1973,56 (7),65-68. Dalle Molle, D. T.; Himmelblau, D. M. Fault Detection in a SingleStage Evaporator via Parameter Estimation Using the Kalman Filter. Ind. Eng. Chem. Res. 1987,26,2482-2489. Frank, P. M. Advanced Fault Detection and Isolation Schemes Using Nonlinear and Robust Observers. Presented at the IFAC 10th Triennial World Congress, Munich, 1987;pp 63-68. Gustafson, R. R.; Slelcher, C. A.; Mckean, W. T.; Finlayson, B. A. Theoritical Model of the Kraft Pulping Process. Znd. Eng. Chem. Process Des. Deo. 1983,22,87-96. Hatton, J. V. Application of Empirical Equations to Kraft Process Control. Tappi 1973,56(€9,108-111. Hatton, J. V.; Keays, J. L. Relationships Between Pulp Yield and Permanganate Number for Kraft Pulps: I. Western Hemlock, White Spruce and Lodgepole Pine. Pulp Pap. Mag. Can. 1970, 71 (11-12),123-132. Hinricks, D. D. The Effect of Kraft Pulping Variables on Delignification. Tappi 1967,50(41,173-175. Jo, J. H.; Bankoff, S. G. Digital Monitoring and Estimation of Polymerisation Reactors. AIChE J. 1976,22,361-368. Jutila, E. A. A. A Survey of Kraft Cooking Control Models. In Instrumentation and Automation in the Paper, Rubber, Plastics and Polymerization Industries; Van Cauwenberghe, A.; Ed.; Pergamon Press: Oxford, 1980. Kuester, J. L.; Mize, J. H. Optimization Techniques with Fortran; McGraw-Hill: New York, 1973;Chapter 9. Macdonald. The Puloina . - of. Wood:McGraw-Hill: New York, 1969; Vol. 1, Chapter 9. Park. S.: Himmelblau. D. M. Fault Detection and Diagnosis via Parameter Estimation in Lumped Dynamic Systems. Ind. Eng.
Ind. Eng. Chem. Res. 1992,31,855-860 Chem. Process Des. Dev. 1983,22 (3), 482-487. Ramirez, W.F. Optimal State and Parameter Identification. An Application to Batch Beer Fermentation. Chem. Eng. Sci. 1987, 42, 2749-2756. Shannon, R. E. Systems Simulation: The Art and Science; Prentice-Hall: Englewood Cliffs, NJ, 1975;p 371. Stephanopouloa,G.; San, K. Y.Studies on On-line Bioreactor Identification: Theory. Biotechnol. Bioeng. 1984,26,1176-1188. Tappi Standards and Suggested Methods; Technical k i a t i o n of the Pulp and Paper Industry: 360 Lexington Avenue, New York, 1962.
855
Venkateswarlu, Ch.; Gangiah, K.; Jatkar, D. D.; Sridharan, S. S. Modelling and Simulation of Kraft Pulping Digester. Indian J. Technol. 1986,24,667-670. Vroom, K.E. The H-factor: A Means of Expressing Cooking Times and Temperatures as a Single Variable. Pulp Pap. Mag. Can. 1957,58,228-231. Wells, C. H. Application of Modern Estimation and Identification Techniques to Chemical Processes. AZChE J. 1971,17,966-973.
Received for review September 9, 1991 Accepted October 30,1991
Derivation of Transfer Function from Relay Feedback Systems Ren-Chiou Chang, Shih-Haur Shen, and Cheng-Ching Yu* Department of Chemical Engineering, National Taiwan Institute of Technology, Taipei, Taiwan 10772, ROC
Luyben proposed a procedure for system identification based on the relay feedback system. This approach can lead to significant errors in the system parameters for transfer functions with time delay. Analytical expressions for the period of oscillation and amplitude ratios are derived. On the basis of these expressions, a procedure is proposed for deriving transfer functions. The results show that improvement can be achieved using the proposed method. 1. Introduction On the basis of the hom-Hagglund (1984) relay feedback system (ATV, autotune variation),Luyben (1987) proposed a procedure for the identification of process transfer functions. The procedure consists of two steps. First, a relay feedback experiment is conducted, and the ultimate gain and ultimate frequency are recorded. Second, this information is fitted to a typical transfer function in chemical process control, e.g., first-, second- or thirdorder plus time delay system. This identification approach, called the original ATV method hereafter, becomes a standard practice in chemical process control, as can be seen in recent textbooks in process control (Seborg et al., 1989, p 301; Luyben, 1990,p 520). The distinct advantages of the ATV method are as follows: (1)it identifies process information around the important frequency, ultimate ; frequency (the frequency where the phase angle is a)(2) it is a closed-loop test, and therefore the process will not drift away from the nominal operating point; and (3) it is a more efficient method than conventional step or pulse testing (requires less time in the experiment). Despite the apparent success in industrial applications, the ATV method can be improved significantly to obtain a more accurate estimate of the system papmeters. The reason is that the estimated ultimate gain K,, and estimated ultimate frequency ij,, derived from the describing function are only an approximation of information at the critical frequency. This approximation can lead to significant errors in the ultimate gain K,,and ultimate frequency a,,(e.g., 5-20% error in K,,)for typical transfer functions in process control systems. The purpose of this work is to derive analytical expressions for the period of oscillation and amplitude ratio for the relay feedback systems such that a better estimation of the process transfer function can be achieved. This paper is organized as follows. Section 2 reviews the ATV method of Luyben (1987). Analytical expressions for the period of oscillation and amplitude ratio are derived for typical process transfer functions (models 1-5 of Luyben (1987)),and an identification procedure based on the new results is presented in section 3. Examples are used to illustrate the accuracy of the proposed procedure, and the
*Towhom all correspondence should be addressed.
robustness with respect to errors in the time delay is investigated in section 4, followed by Conclusions. 2. Review of ATV Method The Astrom-Hiigglund relay feedback system is based on the observation that a feedback system in which the output lags behind the input by -T rad may oscillate with a period P,,.Figure 1illustrates the relay feedback system. A relay of magnitude h is inserted in the feedback loop. Initially, the input u is increased by h. As the output y starta to increase (after a time delay D), the relay decreases the input to a position h below the steady-state valus. Since the phase lag is -T, a limit cycle with a period P,, results (Figure 1). The period of the limit cycle is the ultimate period. Therefore, the estimated ultimate frequency from this relay feedback experiment is iju = 2 4 , , (1) From the Fourier series expansion, the amplitude a can be considered as the result of the primary harmonic of the relay output. Therefore, the ultimate gain can be approximated as (Ogata, 1970; Astrom and Hagglund, 1984) K,,= 4 h / ~ a (2) After the reJay feedback experiment, the estimated ultimate gain ( K J ,estimated ultimate frequency (ij,,), and the observed time delay are used to back-calculate the process transfer function. Typical transfer functions in process control are assumed and parameters are calculated. The transfer functions have the following forms (Luyben, 1987): models 1-3 Kpe-Ds n = 1, 2, 3 G(s) = (3) (7s 1)" model 4 Kpe-Ds G(s) = (4) (71s 1)(Tfi+ 1) model 5 Kpe-Ds G(s) = (5)
+
(71s
+ 1)2(Tfi + 1)
0888-5885/92/2631-O855$O3.OO/O 0 1992 American Chemical Society
