Ind. Eng. Chem. Res. 2003, 42, 4043-4049
4043
Incorporating Prior Knowledge in a Cubic Spline ApproximationsApplication to the Identification of Reaction Kinetic Models Ja´ nos Mada´ r,† Ja´ nos Abonyi,*,† Hans Roubos,‡ and Ferenc Szeifert† Department of Process Engineering, University of Veszprem, P.O. Box 158, H-8201 Veszprem, Hungary, and Control Systems Engineering, Delft University of Technology, P.O. Box 5031, 2600 GA Delft, The Netherlands
Data smoothing and resampling are often necessary in the handling of data obtained from laboratory and industrial experiments. This paper presents a new algorithm for incorporating prior knowledge into the spline smoothing of interrelated multivariate data. Prior knowledge based on visual inspection of the variables and/or knowledge about the assumed balance equations can be transformed into linear equality and inequality constraints on the parameters of the splines. The splines can then be simultaneously identified from the available data by solving one quadratic programming problem. To demonstrate the applicability of this method, two examples are given. In the first example, the proposed approach is applied to the identification of the kinetic parameters of a simulated reaction network, whereas in the second example, data taken from an industrial batch reactor are analyzed. The results show that, when the proposed constrained spline-smoothing algorithm is applied, one achieves not only better fitting to the data points, but also improved performance in the estimation of the kinetic parameters compared to the case in which no prior knowledge is involved. 1. Introduction In many practical situations, laboratory and industrial experiments are expensive and time-consuming, and exact measurements cannot be made. This problem results in small numbers of data points that are often noisy and obtained at irregular time intervals. Hence, data smoothing and resampling are often required to handle such data sets. A good method for this purpose is the cubic spline approximation,1,2 which provides acceptable results and has practical relevance;3 for example, it is widely used for the estimation of reaction kinetic models.4,5 The identification of reaction rates is an important problem, as chemical process models often contain reaction networks, and the kinetic parameters of these models often cannot be identified on the basis of a priori knowledge alone. Therefore, process modeling is generally supported by experiments to identify the kinetic parameters. Many methods have been suggested to obtain reasonable estimates for these rate coefficients from experimental data.6 For example, in the well-known paper of Himmelblau and co-workers, an interesting method for estimating rate constants of complex kinetic models from isothermal batch or plug-flow reactor data was described.7 This method, similarly to other approaches, is applicable when large numbers of data points along the concentration trajectories are available. To avoid the problem of obtaining large collections of data, individual splines are fitted for each measured variable, and the data are resampled and smoothed on the basis of these splines.3-5 Unfortunately, the quality of a spline approximation is also severely affected by the number of * To whom correspondence should be addressed. E-mail:
[email protected]. Internet: www.fmt.vein.hu/softcomp. † University of Veszprem. ‡ Delft University of Technology.
the measured data points and the measurement noise. Hence, there is a need for a new approach that can handle such problems. Recently, combinations of a priori knowledge with data-driven modeling techniques have gained considerable interest in an effort to employ as much knowledge as possible by simultaneously using prior-knowledgeand measurement-based information.8 In general, different modeling paradigms should be used for the efficient utilization of these various sources of information. This means that, if reliable mechanistic knowledge about the process is available, such knowledge can be transformed into a white-box model described by analytical (differential) equations. If the information available is based on human experience and described by linguistic rules and variables, the mechanistic modeling approach is useless, and the application of qualitative9 and rule-based approaches such as fuzzy logic is more appropriate.10 Finally, situations might arise in which the most valuable information comes from input-output data collected during process operation. In this case, the application of black-box models is the best choice. Unfortunately, the actual situation is not clearly one of the previously mentioned approaches. In such cases, the success of the modeling effort hinges upon the ability to combine different types of information in a creative, synergistic way. According to the type and amount of available information, four basic levels of model synthesis can be defined:8 (1) A white-box model is a complete mechanistic model that is constructed from a priori knowledge and physical insight. Because typical white-box models consist of balance equations, the parameters and variables of the model have physical meaning. (2) Blackblock models are designed entirely from data. The parameters of these models are tuned to fit the observed data and do not have any physical meaning. (3) A hybrid or semi-mechanistic model, also called an “incomplete
10.1021/ie0205445 CCC: $25.00 © 2003 American Chemical Society Published on Web 07/29/2003
4044
Ind. Eng. Chem. Res., Vol. 42, No. 17, 2003
white-box model” is constructed from a priori knowledge about the physical background of the process but also includes black-box elements to represent the unknown parts of the system.11 (4) Grey-box models are based on black-box model and incorporate prior knowledge in the selection of the model structure and in the parameter identification procedure. Because the spline approximation is a black-box modeling approach, when prior knowledge is incorporated into the identification of the parameters of the splines, a gray-box modeling scheme results. The most important direction of gray-box modeling applies constraints to the model parameters. For example, Timmons et al.12 applied such a gray-box strategy for adaptive control, and Abonyi et al.13 presented a graybox approach in which the a priori knowledge was translated into inequality constraints on the fuzzy model parameters. The main contribution of this paper is the development of a multivariate gray-box spline-smoothing method. The main idea is to constrain the candidate parameters of the splines. As will be shown, prior knowledge based on the assumed material balance and/ or visual inspection of the data can easily be transformed into linear parameter constraints. The splines can then be identified from the available data by quadratic programming. The proposed approach is applied to a simulated data set and a real-world data set taken from an industrial batch reactor. Splines are first obtained from concentration measurements, and the estimated concentration trajectories are then used in the identification of reaction kinetic parameters. The results show that, when the proposed constrained spline-smoothing algorithm is applied, one achieves not only better fitting to the data points, but also improved performance in the estimation of the kinetic parameters compared to the case in which no prior knowledge is involved The paper is organized as follows: In section 2, the applied spline model is detailed. Section 3 describes how prior knowledge can be implemented as constraints on the parameters of the splines. Section 4 presents the application examples. Conclusions are given in section 5. 2. Standard Cubic Spline Approximation Concentration trajectories involving noisy data and large time intervals are often difficult to approximate by low-order polynomials. However, cubic splines, which are piecewise third-order polynomials, seem to provide a suitable alternative. In cubic splines, the polynomials are defined such that their values and first derivatives are continuous at so-called knots where the individual polynomials are interconnected. When such splines are identified, a continuous function is fitted to the available measured data, x ) [x1, ..., xN]T, given at time instants t ) [t1, ..., tN]T. An efficient way to calculate the coefficients of the cubic spline is given by Horiuchi.4 To formulate this algorithm, let us define a cubic spline for the knot sequence t1 ) k1 < k2 < k3 < ... < kn-1 < kn ) tN. The cubic spline is a sequence of cubic polynomials defined for each interval, [k1 k2], [k2 k3], ..., [kn-1 kn], by the combination of the function values and the first-order derivatives at the knots
S(t) ) s′iai(t) + s′i+1bi(t) + sici(t)+ si+1di(t)
(1)
for ki e t < ki+1, where
si ) S(ki), ai(t) )
s′i )
dS(t) | dx t)ki
(2)
(ki+1 - t)2(t - ki) , hi2 bi(t) ) -
(ki+1 - t)(t - ki)2 hi2
(3a,b)
(ki+1 - t)2(2(t - ki) + hi) , ci(t) ) hi3 di(t) )
(t - ki)2(2(ki+1 - t) + hi) hi3
(3c,d)
with hi ) ki+1 - ki. As can be seen from eq 1, the spline is linear in the parameters
θ ) (s1, s′1, s2, s′2, ..., sn, s′n)T Hence, the θ parameter vector can be determined by minimizing the following quadratic cost function
min Q(θ)
where
θ
Q(θ) )
(
1
N
∑[xj - S(tj)]2 ) ||x - Ψθ||2 Nj)1
where Ψ) c(t1) a(t1) d(t1) b(t1) 0 0 ... 0 c(t2) ... b(t2) 0 0 ... 0 ... ... ... ... c(tp) ... b(tp) 0 0 ... 0 0 0 ... 0 ... c(tN) ... b(tN) 0 0 ... 0
)
(4)
(N × 2n)
(5)
This optimization problem can be solved analytically by the ordinary linear least-squares (LS) method to give
θ ) (ΨΨT)-1ΨTx
(6)
The advantage of the utilized cubic spline is that the integral and derivative of the spline are also linear in the parameters of the spline; that is
dS(t) ) s′ia′i(t) + s′i+1b′i(t) + sic′i(t) + si+1d′i(t) dt
(7)
where
a′i(t) )
-2(ki+1 - t)(t - ki) + (ki+1 - t)2 hi2
(8a)
Ind. Eng. Chem. Res., Vol. 42, No. 17, 2003 4045
b′i(t) ) -
c′i(t) )
2(ki+1 - t)(t - ki) - (t - ki)2 hi2
-2(ki+1 - t)(2(t - ki) + hi) + 2(ki+1 - t)2
d′i(t) )
hi3 2(t - ki)(2(ki+1 - t) + hi) - 2(t - ki)2 hi3
Cj,i
(8b)
(8c)
l
Q(θ˜ ) )
∑ j)1
{
1
Nj
∑
Nj i)1
min Q(θ˜ )
} ∑{ l
[xj,i - Sj(tj,i)]2 + λ1
Aj,iSj(tˆi)]2 + λ2
j)1
N ˆ
1
N ˆ i)1
1
N ˆ
∑
[bj,i N ˆ i)1 dSj(tˆi) 2 dj,i - Cj,i (14) dt
} ∑{ ∑[ l
The main idea of this paper is the simultaneous spline smoothing of multivariate data. In this case, the unknown parameter vector contains the parameters of the individual splines of l variables, θ˜ ) [θ1T, ..., θlT]T, where θj ) [sj,1, s′j,1, sj,2, s′j,2, ..., sj,n, s′j,nj]T. The available measured data are arranged into vectors, x˜ ) [x1T, ..., xlT]T and t˜ ) [t1T, ..., tlT]T, where xj ) [xj,1, ..., xj,Nj]T. (Note that Nj does not have to be identical to Ni.) The θ˜ parameter vector can also be estimated by minimizing the following quadratic cost function
(13)
where Sj is the jth spline function and ˆti is an arbitrary time instant used for the evaluation of the constraints, tˆ ) [tˆ1, ..., ˆtNˆ ]T. With the use of these additional equations, the resulting LS estimation problem is formulated as
(8d)
3. Simultaneous Spline Approximation for Several Compounds
dSj(tˆi) ) dj,i dt
j)1
]}
˜ 1θ˜ ||2 + λ2||y˜ 2 - Ψ ˜ 2θ˜ ||2 Q(θ˜ ) ) ||x˜ - Ψ ˜ θ˜ ||2 + λ1||y˜ 1 - Ψ (15) Generally, several constraints can be defined in the form p
Q(θ˜ ) ) ||x˜ - Ψ ˜ 0θ˜ ||2 +
λi||y˜ i - Ψ ˜ iθ˜ ||2 ∑ i)1
(16)
This LS formulation allows for the incorporation of available prior knowledge into the simultaneous spline fitting procedure. Because the parameters of the splines are the function values and the function derivatives at the knots, prior knowledge can be easily represented by equality and inequality constraints defined on the parameters of the model. It is advantageous to apply linear inequality and equality constraints, that is
where the λ1, λ2, ..., λp regularization parameters determine the weights of the soft constraints. This procedure resembles ridge regression,15 where the optimum of the regularization parameters is found by cross-validation. When more than one data set is available, this cross-validation method can also be followed to design these parameters. However, in this paper, a heuristic approach is followed, where the parameters were set according to the belief of the modeler in the accuracy and importance of the corresponding constraints. Certainly, the soft-constrained method can be combined with hard constraints, for example, by using. the soft constraints to regularize (smooth) the spline while using the hard constraints to incorporate the prior knowledge into the parameter identification.
min Q(θ˜ )
4. Application Examples
θ˜
where l
Q(θ˜ ) )
θ˜
∑ j)1
{∑ 1
}
Nj
Nji)1
[xj,i - Sj(tj,i)]2 ) ||x˜ - Ψ ˜ θ˜ ||2 (9)
Gθ˜ e g and Hθ˜ ) h, (10)
subject to
as the resulting quadratic programming (QP) optimization problem has a global optimum and can be effectively solved by standard numerical packages, such as MATLAB.14 Such prior knowledge can be based on visual inspection of the data (e.g., one might observe that the concentration of a component is monotonically decreasing), or it can come from balance equations such as
Wx ) q
(11)
dx )0 dt
(12)
V
The drawback of this hard-constrained method is that it constrains the values and the derivatives of the splines only at the knot points. To constrain the complete spline, a “soft-constrained” approach has been developed that handles equality constraints formulated as
Aj,iSj(tˆi) ) bj,i
To demonstrate the applicability of the above-described method, two examples are presented. In the first example, the proposed approach is used to identify the kinetic parameters of a simulated reaction network, whereas in the second example, it is applied to the analysis of data taken from an industrial batch reactor. 4.1. Application to the Estimation of Rate Constants for Kinetic Models. In this section, the proposed spline-smoothing approach is used to identify the rate constants of kinetic models from concentration profiles. The assumed reaction network is given by the following equation
dCi dt
n
)
wi,jkjRj ∑ j)1
(17)
where t is the time or a time-like variable, wi,j is the stoichiometric coefficient for component i in reaction j, kj is the rate constant for reaction j, and the factors Rj are functions of the concentrations. Using t1 as the initial time, the preceding equation can be integrated to yield
4046
Ind. Eng. Chem. Res., Vol. 42, No. 17, 2003 Table 1. Data Sets
n
Ci,k - Ci,1 )
wi,jkjR h j,k ∑ j)1
(18)
where
R h j,k )
∫tt Rj dt k
1
(19)
If the Ci,k values are known and the R h j,k values can be estimated, eq 18 becomes a system of linear equations in the rate constants ki. Hence, the rate constants can be estimated by ordinary least-squares methods, provided that the number of independent equations is not less than the number of rate constants. When the concentration trajectories are sampled frequently, eq 19 can be integrated by numerical quadrature using the data points.7 This scheme, however, is no longer available when the number of data points is small, because of a lack of precision. In this case, a spline-smoothing procedure should be followed. For example, in the paper of Tang,16 a natural cubic spline was fitted for each Rj. In this example, a new approach is presented in which the splines of the concentration profiles are simultaneously fitted, and constraints based on the component balance equations are incorporated to ensure the consistency of the parameters of the estimated splines. The numerical example is taken from Himmelblau et al.7 and Tang16 and involves a consideration of the following isothermal reactions k1
k3
2
4
} B {\ }C A {\ k k
number of data points
relative noise level (%)
independent noise
1 2 3 4 5 6 7 8 9 10
31 31 31 31 31 12 12 12 12 12
0 5 5 10 10 0 5 5 10 10
yes no yes no yes no yes no
(20) Figure 1. Spline approximations for 12 noisy measurements (data set 9): O, measured concentration of A; 0, measured concentration of B; 4, measured concentration of C; solid lines, real, but unknown concentration trajectories; dashed lines, spline estimations with model 3; dotted lines, spline approximations with model 4 (proposed model).
described by the differential equations
dCA ) -k1CA + k2CB dt dCB ) k1CA - k2CB - k3CB + k4CC dt
data set
(21)
dCC ) k3CB - k4CC dt The initial concentrations at t0 ) 0 are CA ) 1, CB ) 0, and CC ) 0, and the real parameters are k1 ) 1, k2 ) 0.5, k3 ) 10, and k4 ) 5. The aim of the experiment is to identify the kinetic parameters (k1, ..., k4) from the available concentration measurements. According to the quality and quantity of the available data, 10 data sets were used to demonstrate the performances of different methods (see Table 1). These data sets were defined from combinations of different (a) numbers of data points (the data points are taken at 12 or 31 time instants), (b) levels of normally distributed noise added to the measured data (0, 5, or 10%), and (c) types of noise (independent or correlated). Four different approaches were followed for the estimation of the kinetic parameters from these data sets. Model 1. Nonlinear, Direct Search Method. The Nelder-Mead simplex method is used for the identification of the kinetic parameters. The concentration trajectories are calculated by solving the differential equations of the model (eqs 21). This method is very effective, but requires a nonlinear optimization algorithm, which is rather slow and can get stuck in local minima.
Model 2. Raw-Data-Based Estimation. In this case, the raw data points are integrated numerically, and the estimation problem is solved by the leastsquares method. Model 3. Spline-Based Estimation. The standard cubic spline approximation is used, and the smoothed concentration trajectories represented by the splines are analytically integrated to estimate the R h j,k parameters needed for the least-squares method. Model 4. Grey-Box Spline-Based Estimation. This approach is similar to model 3, but it utilizes the following prior-knowledge-based constraints
CA + CB + CC ) 1 dCA dCB dCC + + )0 dt dt dt
(22)
dCB dCC dCA < 0, > 0, >0 dt dt dt where the equality constraints represent the component balance equations and the inequality constraints describe the monotonic changes in the concentrations (see Figure 1). In the cubic spline approximation, the choice of knot sequence plays an important role. The knots can be determined automatically by segmenting the time series defined by the concentration profiles;17 heuristically by visually inspecting the variables, looking for significant
Ind. Eng. Chem. Res., Vol. 42, No. 17, 2003 4047 Table 2. Estimation of Kinetic Parameters (P1)a model
a
number of knots
data set 5
1
2
3
4
6
7
8
9
10
1
0.0023
0.18
0.088
0.37
0.18
0.060
0.30
0.16
0.65
0.31
2
0.024
0.38
0.32
0.75
0.64
0.17
0.42
0.36
0.76
0.64
3
4 3 2
0.035 0.12 0.23
0.34 0.33 0.34
0.29 0.26 0.28
0.68 0.62 0.57
0.58 0.46 0.38
0.064 0.15 0.33
0.40 0.41 0.53
0.34 0.36 0.42
0.78 0.78 0.94
0.67 0.64 0.59
4
4 3 2
0.029 0.13 0.42
0.19 0.24 0.49
0.11 0.16 0.42
0.39 0.44 0.72
0.23 0.25 0.46
0.064 0.15 0.45
0.31 0.34 0.67
0.20 0.24 0.48
0.67 0.68 1.43
0.41 0.41 0.62
Means of 100 independent runs.
Table 3. Estimation of Concentration Trajectories (P2)a model
number of knots
1
a
data set 5
1
2
3
4
6
7
8
9
10
1.5 × 10-7
0.0041
0.0011
0.0016
0.0043
0.003
0.011
0.0033
0.043
0.012
3
4 3 2
0.0003 0.0016 0.0036
0.029 0.023 0.018
0.030 0.022 0.018
0.12 0.086 0.062
0.12 0.082 0.061
0.0005 0.0028 0.0065
0.077 0.056 0.042
0.076 0.054 0.038
0.31 0.22 0.15
0.30 0.21 0.13
4
4 3 2
0.0004 0.0035 0.0117
0.016 0.013 0.018
0.0075 0.0066 0.013
0.059 0.040 0.035
0.0005 0.0035 0.012
0.0005 0.0035 0.012
0.043 0.031 0.027
0.019 0.013 0.016
0.016 0.11 0.069
0.070 0.039 0.028
Means of 100 independent runs.
breaks in the concentration trajectories; or intuitively by using some prior knowledge about the system to be modeled (see the second example in section 4.2). In this example, three equidistant knot sequences with two, three, and four knots were investigated and compared. Splines with more than four knots were not employed because the small number of data points, 12 (data sets 6-10), meant that too few data points would fall in the knot intervals, which would make the estimation problem ill-conditioned. The performances of the models were measured according to the following normalized cost function
P1 )
x∑( ) 4
ki - kˆ i
i)1
ki
2
(23)
where kˆ i represents the estimated kinetic parameters and ki represents the real parameters. Models 1, 3, and 4 can also be evaluated in terms of another cost function, which measures the error between the “real” and the estimated concentration trajectories 3
{∫0 ∑ i)1
P2 ) 100
tN
estimated [creal (t)]2 dt} i (t) - ci
(24)
Tables 2 and 3 show that, when the rate coefficients are estimated from the R h j,k values by the least-squares method, usually more accurate estimates are obtained when splines are used (model 2 vs models 3 and 4). When the constrained spline approximation is applied (model 4), not only are the estimated rate coefficients much more accurate, but also better approximations of the concentration trajectories are obtained (model 3 vs model 4). The grey-box spline estimation method (model 4) always estimates the rata coefficients better than the simple LS method (model 2). The nonlinear method (model 1) offers the best results, but it is much slower than the other methods (see Table 4). Furthermore, model 1 is based on the differential equations of the
Table 4. Computational Effort in the Identification of the Modelsa
a
model
time (s)
1 2 3 4
12 0.0069 0.65 0.67
One typical run on AMD Duron 800-MHz computer.
process. Hence, in contrast to the spline approximation, model 1 requires full knowledge of the differential equations. Thus, if only an approximation of the concentration trajectories is the purpose of the modeler, the application of the proposed spline approach would be recommended. Furthermore, for the estimation of the kinetic parameters, the spline-based approach can be used to provide good initial conditions for the identification of model 1, as was suggested in ref 16. When the measurement noise was correlated (data sets 3, 5, 8, and 10), usually better estimations were obtained compared to the uncorrelated case (data sets 2, 4, 7, and 9). This suggests that the concentration trajectories with correlated noise are less in conflict with each other and with the measured data, so that the constraints can more easily handle such measurement errors compared to the totally random uncorrelated case. Generally, an increase of the number of knots improves the accuracy of the model. However, as with other estimation problems, overfitting can occur, so the number of parameters has an optimal value. As Tables 2 and 3 show, this optimum is a function of the noise level. 4.2. Industrial Batch Reactor. The second example is taken from one of our industrial partners. In the studied industrial batch reactor, only the concentrations of the three main components are measured (see Figure 2). As is shown in this example, the proposed gray-box spline-smoothing approach can be effectively used in such typical real-world situations. Two approaches were used and compared: the standard cubic spline approximation and the simultaneous spline approximation with hard and soft constraints.
4048
Ind. Eng. Chem. Res., Vol. 42, No. 17, 2003
The component balance can be formulated as
CA(t) + CB(t) + CC(t) + CA*(t) + CB*(t) + CC*(t) ) C0 (30) where C0 ) CA0 + CB0 + CC0. Although the concentrations of A*, B*, and C* are not measured, with the use of the assumed kinetic models
dCA*(t) ) k1CA(t) dt dCB*(t) ) k2CB(t) dt
(31)
dCC*(t) ) k3CC(t) dt the material balance can be used in practice in the following form Figure 2. Spline approximation (second example) without contraints: O, measured concentration of A; 0, measured concentration of B; 4, measured concentration of C; solid lines, splines; dashed line, mass balance (A + B + C + A* + B* + C*); dotted lines, calculated concentrations of byproducts (A*, B*, C*).
Because the chemistry of this technology is very complex and not fully known, the reaction network is modeled by the following scheme
A+BfC k1
A 98 A* k2
B 98 B*
(25)
k3
C 98 C* where the first reaction represents product formation and the last three reactions represent the decay of three measured components (byproduct formation). Because the technology is confidential, the components are denoted as A, B, and C. As can be seen in Figure 2, the rates of the reactions are relatively small until 60 min of operation. At about this time, a catalyst was added into the reactor that makes the product formation much faster. This prior knowledge is used in the selection of the knot sequence. Several alternatives were tested with larger and smaller numbers of intervals and different knot sequences, and finally, on the basis of the recipe of the technology and an analysis of the shapes of the trajectories, the knot sequence was set as [0, 57.3, 80.3, 250]. From visual inspection of the data (see Figure 2), the following constraints can be defined
CA(0) ) CA0, CB(0) ) CB0, CC(0) ) CC0
(26)
∫0tCA(τ) dτ + t t k2∫0 CB(τ) dτ + k3∫0 CC(τ) dτ
C0 ) CA(t) + CB(t) + CC(t) + k1
Equation 32 is linear in the parameters k1, k2, and k3 and linear in the θ˜ parameters of the splines. Hence, as the CA(t), CB(t), and CC(t) concentration trajectories are represented by the S1, S2, and S3 splines, the product of these parameters appears in the last three terms of eq 32. Hence, this material balance represents a bilinear constraint defined on the unknown parameters. Such bilinear optimization problems are often solved by alternating optimization (see ref 18 for a practical application in system identification). This approach results in the following algorithm: Step 1. Hard-Constrained Spline Identification. The CA(t), CB(t), and CC(t) concentration trajectories are represented by the S1, S2, and S3 spline functions obtained by simultaneous spline approximation using hard constraints, eqs 26-29. Step 2. Identification of the Parameters k1, k2, and k3. Using eq 32, the three kinetic parameters can be identified by the ordinary least-squares algorithm
min Q(k1,...,k3) where
k1,..,k3
N ˆ
Q(k1,k2,k3) )
∑ i)1
3
{C0 - [
∑ j)1
θ˜
(27)
dCC(t) g 0 if t e 114 dt
(28)
dCA(t) dCB(t) dCC(t) + + e0 dt dt dt
(29)
kj∫0 ∑ j)1
ˆti
Sj(τ) dτ]}2
Step 3. Soft-Constrained Spline Identification. In this step, we incorporate the material balance in eq 32 into the spline identification as a soft constraint, so the new quadratic objective function is
minQ(θ˜ ) where Q(θ˜ ) ) dCB(t) dCA(t) e 0, e0 dt dt
3
Sj(tˆi) +
(33)
3
and
(32)
λ
1
N ˆ
3
{
1
Nj
∑ ∑[xj,i - Sj(ti)]2 j)1 N i)1 j 3
Sj(tˆi) + ∑kj∫0 ∑{C0 - [∑ N ˆ i)1 j)1 j)1
ˆti
}
+
Sj(τ) dτ]}2 (34)
Because Sj(t) and ∫t0Sj(τ) dτ are linear functions of the θ˜ parameters and minθ˜ Q(θ˜ ) is subjected to the constraints in eqs 26-29, the optimization problem can be solved by quadratic programming. The results of this
Ind. Eng. Chem. Res., Vol. 42, No. 17, 2003 4049
step are the S1, S2, and S3 splines related to the CA(t), CB(t), and CC(t) concentration trajectories. Step 4. Iteration. Repeat from step 2 until the estimated trajectories converge. Before the application of the soft constraint, the tˆ ) [tˆ1, ..., ˆtNˆ ] time instants used for the evaluation of the soft constraints and the λ parameter have to be selected. In this example, the soft constraints were evaluated in every minute of the operation, and the λ parameter was set according to a rough sensitivity analysis. An increase in λ increased the accuracy of the mass balance. When λ was greater than 15, this increase did not cause further significant improvement, so this parameter was chosen as 15 (see Figure 3).
tions) and the visual inspection of the data (e.g., monotonicity). The proposed approach has been applied to the estimations of concentration trajectories and the parameters of chemical reaction kinetic models. The presented examples showed that, when the proposed identification approach is applied, one obtains not only better approximations of the concentration trajectories, but also more accurate estimates of the rate coefficients. Acknowledgment The authors acknowledge the support of the Cooperative Research Center (VIKKK) (KKK-I-7), the Hungarian Ministry of Education (FKFP-0073/2001), and the Hungarian Science Foundation (OTKA TO37600). J.A. is grateful for the financial support of the Janos Bolyai Research Fellowship of the Hungarian Academy of Science. Literature Cited
Figure 3. Spline approximation (second example) with constraints: O, measured concentration of A; 0, measured concentration of B; 4, measured concentration of C; solid lines, splines; dashed line, mass balance (A + B + C + A* + B* + C*); dotted lines, calculated concentrations of byproducts (A*, B*, C*).
Table 5 shows some comparisons of the numerical results of the two methods. The proposed constrained method fulfils the mass balance very well, in contrast to the standard spline approximation method. Certainly, the spline that was calculated by the standard method fits the measured data better, but the difference between the two methods is quite small. Hence, these numerical results support the conclusion that the proposed method is more accurate. Table 5. Numerical Results of Spline Methods in the Second Example
mean absolute difference between estimated and measured concentrations (g/dm3) mean absolute error in the mass balance (g/dm3) absolute error in the balance at the end of the experiment (g/dm3)
without constraints
with constraints
5.29
5.99
6.36
0.18
10.50
0.08
5. Conclusions In this paper, a constrained spline approximation method has been developed to incorporate prior knowledge into the preprocessing of interrelated, noisy, and rarely sampled process data. The constraints are based on the available prior information (e.g., balance equa-
(1) Ahlberg, J. H.; Nilson, N. E.; Walsh, J. L. The Theory of Splines and Their Application; Academic Press: New York, 1967. (2) Shikin, E. V.; Plis, E. V. Handbook on Splines for the User; CRC Press: Boca Raton, FL, 1995 (3) Kolla´r-Hunek, K.; La´ng-La´zi, M.; Herrmann, N.; Miklo´s, D.; Kova´cs, I. Application of special numerical approximation in thermodynamics. Hung. J. Ind. Chem. 1998, 26, 269-274. (4) Horiuchi, J.; Kishimoto, M.; Kamasawa, M.; Miyakawa, H. Data Based Simulation of Batch and Fed-Batch Cultures for R-Amylase Production Using a Culture Data Base and a Statistical Procedure. J. Ferment. Bioeng. 1993, 76 (4), 326-332. (5) Maria, G.; Muntean, O. Model Reduction and Kinetic Parameters Identification for the Methanol Conversion to Olefins. Chem. Eng. Sci. 1987, 42 (6), 1451-1460. (6) Schaich, D.; Becker, R.; King, R. Qualitative Modeling for Automatic Identification of Mathematical Models of Chemical Reaction Systems. Control Eng. Pract. 2001, 9, 1373-1381. (7) Himmelblau, D. M.; Jones, C. R.; Sichoff, K. B. Determination of Rate Constants for Kinetics Models. Ind. Eng. Chem. Fundam. 1967, 6 (4), 539-543. (8) Sjoberg, J.; Zhang, Q.; Ljung, L.; Benvebiste, A.; Delyon, B.; Glornnec, P. Y.; Hjalmarsson, H.; Judutsky, A. Nonlinear blackbox modelling in system identification: A unified overview. Automatica 1995, 31 (12), 1691-1724. (9) Schaich D.; Becker, R.; King, R. Qualitative modelling for automatic identification of mathematical models of chemical reaction systems. Control Eng. Pract. 2001, 9, 1373-1381. (10) Abonyi, J. Fuzzy Model Identification for Control; Birkhauser: Boston, 2003. (11) Psichogios, D. C.; Ungar, L. H. A Hybrid Neural NetworkFirst Principles Approach to Process Modeling. AIChE J. 1992, 38 (10), 1499-1511. (12) Timmons, W. D.; Chizeck, H. J.; Casas, F.; Chankong, V.; Katona, P. G. Parameter Constrained Adaptive Control. Ind. Eng. Chem. Res. 1997, 36, 4894-4905. (13) Abonyi, J.; Babuska, R.; Verbruggen, H. B.; Szeifert, F. Incorporating Prior Knowledge in Fuzzy Model Identification. Int. J. Syst. Sci. 2000, 31 (5), 657-667. (14) MATLAB Optimization Toolbox; MathWorks, Inc.: Natick, MA, 2002. (15) Marquardt, D. W.; Snee, R. D. Ridge regression in practice. Am. Stat. 1975, 29 (1), 3-20. (16) Tang, Y. P. On the Estimation of Rate Constants for Complex Kinetics Models. Ind. Eng. Chem. Fundam. 1971, 10 (2), 321-322. (17) Chatfield, C. The Analysis of Time Series: An Introduction; Chapman & Hall: New York, 1996. (18) Abonyi, J.; Babuska, R.; Botto, M. A.; Szeifert, F.; Nagy, L. Identification and Control of Nonlinear Systems Using Fuzzy Hammerstein Models. Ind. Eng. Chem. Res. 2000, 39, 4302-4314.
Received for review July 22, 2002 Revised manuscript received April 18, 2003 Accepted June 30, 2003 IE0205445
