Znd. Eng. Chem. Res. 1995,34, 567-574
567
Design and Experimental Evaluation of a State Estimator for a Crystallization Process Rob A. Eekt and Sjoerd Dijkstra" Mechanical Engineering Systems and Control Group, Delft University of Technology, Mekelweg 2, 2628 CD Delft, The Netherlands
The development of a state estimator for a crystallization process, equipped with a sensor for the crystal size distribution (CSD), is discussed. The estimator is designed on the basis of a nonlinear distributed parameter model, which describes both the dynamics of the crystal size distribution and the output of the sensor. The process model is lumped to a set of nonlinear first-order differential equations, which are linearized leading to a high-order state-space model. With this model, a constant error feedback gain is designed by using Kalman filter theory. Estimates are used for the measurement and process noise covariance matrices. The designed gain is implemented in the nonlinear lumped process model. The resulting state estimator is evaluated with raw data from a pilot crystallizer, equipped with a n on-line CSD sensor. For different sets of output data, the designed estimator is able to track the process output signal trend, while reconstructing the CSD, the supersaturation level, and a set of related variables. Sufficient robustness is demonstrated for sensor failure, initial state errors, and process disturbances.
1. Introduction The dynamic behavior of continuous industrial crystallizers is often characterized by badly dampened or even oscillatory behavior in time (Eek et al., 1995). The performance of a suspension crystallizer is mainly determined by the size, the purity, and the morphology of the produced crystals. Although crystal purity and morphology are relatively new and active areas of present academic and industrial research, the size of crystals expressed by the crystal size distribution (CSD) is often considered as the most important process parameter to be optimized. Subsequently, the supersaturation of the solution, which enables crystal growth, is seen as an important parameter, since it influences crystal nucleation and growth kinetics, which strongly determine the crystal size, purity, and morphology of particles. In industrial practice, both undesired dynamics of the CSD and the supersaturation limit the marketability and transportability of crystals and put a constraint on the performance of the downstream crystal-handling processes, like thickeners, centrifuges, and dryers. Therefore, a demand exists for control systems that increase the overall performance of crystallizers and downstream processes. A basic requirement for the implementation of a control system is the availability of accurate and reliable information on the relevant process dynamics. For on-line measurement of both the CSD and the supersaturation in crystallizers, different sensors have been developed which, however, still have severe shortcomings. On-line measurement of supersaturation is usually extremely difficult, as the relative level of supersaturation is often not more than 1%. CSD measurement is difficult as most CSD sensors are indirect; i.e., reconstruction schemes are required t o obtain a proper estimate of the CSD based on the raw output of the sensor. Besides, robust implementation
* Author to whom all correspondence should be addressed. E-mail:
[email protected]. Current address: Bayer AG, ZF-TST, Building E41,D-51368 Leverkusen, Germany. 0888-5885/95/2634-0567$09.00/0
of on-line CSD sensors is difficult (Eek et al., 1993). Solutions that are often applied to the CSD estimation problem are direct inversion algorithms which estimate the CSD based on the estimated inverse of the sensor model. This inversion procedure is, however, ill conditioned and as a consequence sensitive t o measurement noise. In addition, the resolution of the estimated state is limited and the inversion algorithms do not take into account the evolution in time of the measured signal trend, governed by the process dynamics. Another approach to estimate the CSD is to use dynamic process information in combination with the sensor output. This can be achieved with state estimators and observes, which were introduced in the 1960s by Kalman and Luenberger. State estimators and observers are based on a process model and estimate the states of the process from available (limited) input and output measurements. Basically, a state estimator is a dynamic system, which runs on a computer and operates sequentially on the raw measured sensor data. An optimally designed estimator will give the best compromise between sensitivity to measurement noise and rapid convergence to the real process outputs. In contrast with state estimators, observers are mostly designed for deterministic systems with insignificant process or sensor disturbances. The purpose of this article is t o present a practical state estimator for a crystallization process, based on a first principles process model. This estimator should primarily reconstruct the CSD and the superaturation level for control purposes. Second, the output of the estimator should serve as a process monitor to facilitate the operator task in making decisions with respect t o manual control actions. Most observer applications reported in the literature deal with space systems or navigation systems. Unfortunately, the number of applications in the chemical process industry is still limited due to severe modeling limitations and sometimes to the absence of accurate and robust on-line measurement systems. Some successful implementations in polymerization reactors are reported by Dimitratos et al. (1989), and in monitoring for hazardous conditions by Gilles and Schuler (1981). 0 1995 American Chemical Society
568 Ind. Eng. Chem. Res., Vol. 34, No. 2,1995
I
Qv
vcfj
C
T P
fines (Tr), and the total heat input (Pt,t), at prescribed set-point values. Effective actuators for the process dynamics are the fines removal rate, the total heat input to the system defined by Pbt = Pi, P,, and the product removal rate. CSD measurements are performed at the sample location (SL), where the CSD of unclassified removed product crystals is measured with a Malvern particle sizer. The CSD is characterized by the population density function n(x,t), defined by
+
hN(x,t)
n(x,t)a lim Ax-+
Ax
(1)
with N(x,t)the cumulative number function describing the number of crystals per unit volume, with a size equal or lower than x . I
1
Qi
Figure 1. Schematic drawing of an evaporative draR tube baffled crystallizer, equipped with a fines removal system.
Applications of state estimation techniques in the field of crystallization are reported by Cooper and Clough (1986), Hashemi and Epstein (1982), Myerson et al. (19871, and Tsuruoka and Randolph (1987). Most of these applications, however, lack robust on-line process sensors and are often based on simple (MSMPR)process models. In this paper an experimental 970 L pilot crystallizer equipped with a robust CSD sensor was used to generate process data. The CSD measurement system comprises a continuous slurry diluter and a sensor that is based on forward light scattering. The estimator is based on a lumped version of a nonlinear distributed parameter model for the CSD dynamics. This lumped model is chosen instead of the distributed model as it strongly facilitates the filter design and implementation (Anderson and Moore, 1979). A detailed description of the process model and the estimation of unknown model parameters is given by Eek et al. (1994). Step response analysis, based on a linearized version of the model around its steady state, indicates that the crystallizer can be considered as a sufficiently linear process. For this reason we applied linear Kalman filter theory to estimate a constant Kalman gain matrix. The resulting filter is tuned and evaluated on experimental data obtained from the pilot crystallizer. Several practical robustness features, such as sensitivity of the estimator for sensor failure and process disturbances, are outlined. The approach we follow is applicable to other crystallization processes, equipped with different sensors, provided that sufficiently accurate models for their behavior are available. 2. Process Description
The crystallizer studied here is an evaporative isothermal continuous draft tube baffled crystallizer equipped with an annular zone around the crystallizer vessel in which classified removal of the fine particles is established. A n elementary process flow sheet with the process inputs, outputs, and some constant controlled internal variables is depicted in Figure 1. PI control loops are present to maintain the slurry temperature (T), the pressure in the crystallizer (Pr), the crystallizer volume (VJ, the fines (Qf) and product (Qp) removal flows, the temperature of the dissolved
3. Process Modeling A widely accepted framework for the modeling of crystallization processes is a population balance model in which the CSD dynamics are described by a firstorder hyperbolic partial differential equation (Randolph, 1988): -an(x,t) =-
at
aGe(x,t)n(x,t)ax
QW ,
+ H&,t)yXt)QLt) n(x,t) (2) vc
with a boundary condition n(x = xo, t ) and an initial condition n,(x,t = to), which often describes an initial population of seed crystals. The function hf describes the probability that small particles are removed from the crystallizer volume. This function is described by an S-shaped curve. The factor yf is a mass separation efficiency factor. The crystal growth rate Ge(x,t)is expressed as the product of a kinetic growth rate and a length-dependent growth rate: Ge(x,t)= G&)G,(x). The kinetic growth rate is related to the level of supersaturation of dissolved crystal material in the clear liquor:
Gk(t)= p@C(t)P'
(3)
The power p7 is often taken equal to 1 (Randolph and Larson, 1988). The supersaturation is calculated from a nonlinear ordinary differential equation, which can be derived from both mass and heat balances: (4) With u ( t ) the process input vector, which is given by u(t) = (QJt), &At),Pbt(t)},and the process state z , given by z(x,t) = (n(x,t),AC(t)}.The length-dependent growth rate function is given by G,(x)
1-
x%,"" x,""(P
+x p ) +x p )
(5)
which is an S-shaped curve, with xa as the point of inflection and xe as the maximum simulated crystal size, both parametrized by the unknown parameters pg and p10, respectively. The sharpness of this curve is parametrized by the constant factor pa.
Ind. Eng. Chem. Res., Vol. 34, NO. 2, 1995 569 product particle flow He-Ne laser U
c
expander
.93
-
.9 -
1
,
--_._..*-,
,
,
' I
Figure 2. Principle of forward light scattering.
The crystal birth rate B is assumed to be dependent on the supersaturation level and the presence of other crystals:
B(t)= p , [ sP4 m n(x,t>xP5dxlp' AC(t)Pz
(6)
with p1, ..., p5 empirical constants. The factor p3 describes the combined effect of stirrer speed, pump rotation rate, slurry temperature, crystallizer geometry, etc. The crystal birth rate together with the crystal growth rate determine the boundary condition n(xo,t) (Randolph and Larson, 1988): n(x0,t) = B(t)/G(t). The empirical relations contain a set of unknown parameters 8 = ( P I , pz, ..., plo}, which are estimated @omprocess data resulting in optimal parameter value 8. Additional details and the results of the parameter estimation procedure are given in Eek et al. (1995). The Sensor Model. The crystallizer model is extended with a sensor model to describe the output of the particle sizer as a function of the CSD that passes through an optical cell. The working principle is outlined in Figure 2. The outputs of this sensor are 31 values for the intensities of scattered laser light, recorded on the 31 concentrical detector rings (Figure 2). Measurement of particle size distributions with lightscattering techniques is a well-known and often applied technique (Boxman, 1990). It has the advantage that it can cover a relatively wide range of crystal sizes and that it is a fast technique with a good detectibility of error sources. Fraunhofer theory is used t o develop a model for the scattering of the incident laser beam by the differently sized crystals that flow through the optical cell. The crystals are assumed t o have a rectangular shape with an aspect ratio a. The equations needed to develop the scatter model are obtained from Hecht (1987). The discretized version of this model is written as
y=Hn
(7)
with H the sensor model matrix, y the sensor output vector, which contains 31 values for the scattered laser light, and n a discrete population density function. Inversion of the model, given by eq 7, can be used to calculate the size distribution from an average intensity vector y. Boxman (1992) developed an algorithm that estimates the CSD from the weighted inverse of the sensor model:
a, = ( H ~ V H ~ H ~ VR, Y>, o
(8)
The nonnegativity constraint is added t o avoid physically impossible solutions. The weighting matrix V is the measurement noise covariance matrix that corresponds to the averaged light intensity vector y. Both the average intensity vector y and the matrix V are calculated, at each sample instant, from a batch of 1000 sweeps which are read from the detector in a short period of time (approximately 15 s). The resulting CSD, R,(x,t), denotes a discrete crystal volume based size
v .7t 0
'
1
2
'
3
4
'
time
'
5
6
'
7
'
8
9
'
1
'
- -' - _ _ 4_
0
1
1
[hours]
Figure 3. Comparing of linear (dashed) and nonlinear responses of three states: n(2), n(40),and n(80),on a step in the product removal rate.
distribution on a grid with logarithmically equidistant size classes. From this distribution characterizing quantities such as a median crystal size or a spread are calculated. We will compare the outcome of this direct inversion approach with the results obtained from the state estimator.
4. Model Lumping and Linearization Lumping of the continuous population density n onto a set of discrete crystal size classes is applied to transform the partial differential equation to a set of nonlinear ordinary differential equations. A secondorder finite difference scheme is applied to approximate the first term in the right-hand side of the population balance (eq 2) as described by De Wolf et al. (1990).Time integration of the lumped population balance is performed with the forward Euler method with truncation error correction. The model equations are linearized analytically by applying a first-order Taylor expansion arround the steady-state solution of the model. The steady-state solution is calculated by putting the righthand side of eqs 2 and 4 to zero. As the resulting equations do not represent an explicit solution for the stationary state, an iteration scheme is applied. The linearized equations constitute a high-order linear state-space model, given by
A2 = Fhz + GAu,
z(t = to), Ay = H h z
(9)
where A denotes a small perturbation around the stationary solution. For ease of notation, the A symbol is omitted from the equations presented below. To improve the numerical condition, an equivalence transformation is applied, by transforming the process state with a diagonal transformation matrix, which has the inverse of the stationary process state on its diagonal. Lumping of the continuous population density n onto 99 discrete equidistant crystal sizes, in the range of [O-15001 pm, gives a sufficiently accurate approximation of the dynamics to be described. Consequently, the transformed linear state-space model, which is used for the filter design, has 100 states, e.g., 99 states for the discretized CSD n and 1for the supersaturation AC. In Figure 3, the step response of the linear model on a positive step of 10%in the product removal rate, which
570 Ind. Eng. Chem. Res., Vol. 34,No. 2, 1995 process
Z .3
Y
L W .
3 0 . 20.
A
Figure 5. Experimental time response of scaled light intensity ( 2 ) on first 20 detector rings ( x ) over 30 h in time Cy) (fines removal rate 1.0 Us). Figure 4. Schematic drawing of a continuous Kalman filter.
is considered as a process input, was compared with the positive and inverse negative nonlinear model response. For a model with negligible nonlinearities, the positive response should equal the inverse negative response. Disturbing the product removal rate is considered as a worst case, because the responses to other disturbances, e.g., the fines flow and the heat input, were found to be quasi-linear. Therefore, it is assumed that the population balance model has only mild nonlinearities, which can be approximated by linear models, sufficiently accurate.
For linear time-invariant state-space models, the solution for the optimal filter gain is given by the solution of an algebraic Ricatti equation (Anderson and Moore, 1979). When it is assumed that the model, described by eq 9, is disturbed with uncorrelated random process and sensor noise inputs, this model can be written as (A symbol is omitted)
&(t) = F d t ) + Gu + w(t), y ( t ) = Hz(t) + v(t) (12) The corresponding covariance matrices for the sensor noise u(t) and the process noise w(t) are given by
5. Process Observability The state-space model {F,G, H) enables the estimation of observability properties. A n important tool to study observability properties of a process is the observability Gramian Q, which is defined as (Moore, 1981)
Q=
eRHTHeRdt
An interpretation of the observability Gramian can be given in terms of energy: the amount of energy measured at the output of the process from t = 0 to t = -, due to an initial condition xo at t = 0, with zero inputs, is
If Q is singular, an unobservable initial condition will exist that will not contribute to the output signal energy. In addition, Q will provide information about well or weakly observable states.
6. Kalman Filter Design
With Kalman filtering theory, a state estimator can be devised that takes into account the stochastic nature of model and process information in a systematic way. The working principle of a Kalman filter is depicted in Figure 4. The error E (also denoted as the innovation process) between the actual process output y and the model output 9 is calculated and amplified with the so called Kalman gain Kf.This results in a correction signal that is used as an additional input to force the model in a way that the estimated state converges to the real process state.
E(vTv) = V
(13)
E(wTw)= W
(14)
Based on this model, the filter gain matrix Kfis given bY
where Pfsatisfies an algebraic Riccati equation:
+
PfFT FPf- PpTTIHPf+ W = 0
(16)
This solution indicates that an accurate process and sensor model, and good estimates for the covariance matrices for the different noise sources, play a decisive role in the filter design procedure. This can be seen as a disadvantage as it is ofien difficult or even impossible to obtain a sufficiently accurate model representation of the process and the corresponding noise sources. However, in many practical situations, it suffices to treat the noise covariance matrices as free design parameters that are adapted to achieve a satisfactory filter response based on an imperfect model. A designed . filter can be denoted as optimal when the error E is white with a zero mean.
7. Results and Discussion
Experimental Data. The results of two free-run experiments (described by Eek et al. (1995)), measured on the product output (SL in Figure 1)of the pilot plant a t different process conditions, were used to optimize and evaluate the designed state estimator. Figures 5 and 6 show the raw light intensity patterns measured on the inner 20 detector rings of the Malvern instrument, over 30 h for a fines removal rate of 1.0 and 3.4 Us,respectively. The measured signal trends are
Ind. Eng. Chem. Res., Vol. 34, No. 2, 1995 571
.
Z .4
r-
.9
-
.8
-
n 0 J
-0
-7
-
L
.6
-
0,
-0J
.5
-
C
X
n
-0
.4
-
Y
-3 -
e
A
30.20.
I:
X
0
I
Figure 6. Experimental time response of scaled light intensity ( 2 ) on first 20 detector rings ( x ) over 30 h in time (y) (finesremoval rate 3.4 Us).
.2 -
I - 0
loo0
I A
.!E
u
200 1000
I.'.'' -
'
-1
- ' ~ ~ . . ' . . . " . ~ . ' ' . ' ~ ' ' '
.
'
I
0 10 X
0
5
20
40
60
80
100
Figure 8. Singular values of the observability Gramian.
10
15
20
25
30
35
40
time [hours]
Figure 7. Experimental time response of X50 (fines removal rate 1.0 and 2.2 Us,respectively).
corrected for outliers by linear interpolation. In Figure 7, the corresponding experimental time response is given for the median crystal size x50, which was obtained from volume-based distributions calculated by direct inversion (eq 8). In this plot all outliers caused by tube blockage or other process disturbances are present. As can be seen, large fines flows induce limit cycles. For this behavior the validity of our approach may be questioned, since it is based on the assumption that a stationary linear model is sufficiently descriptive for the process dynamic behavior arround a stationary state. In the first instance, a more suitable approach for this case is offeredby extended Kalman filtering (Anderson and Moore 1979), which is a straightforward extension of the linear stationary Kalman filter. An extended Kalman filter solves a set of recursive filter equations a t discrete time instances along a state trajectory that is followed by the process. The model matrices used in this approach can be obtained by repeatedly linearization of the nonlinear model at discrete time instants along the followed state trajectory. The resulting Kalman feedback gain will be time varying, however; for small deviations from the stationary state, the filter gain will approach the stationary filter gain. The extended Kalman filter approach was tested on the process data; however, it did not lead to a significant improvement of the solution; it merely resulted in an excessive computational burden. On the basis of this,
and the results presented below, it is concluded that the use of a stationary feedback gain derived from the solution of a stationary Ricatti equation gives satisfactory and applicable results. ObservabilityAnalysis. Inspection of the singular values of the observability Gramians, depicted in Figure 8, reveal that only a limited number of modes (
