Parameter Estimation of Unsteady-State Distributed Models in the Laplace Domain John A. Williams’ and Robert J. Adler Chemical Engineering Science Division, Case Western Reserve University, Cleveland, Ohio
William J. Zolner 111 Department of Chemical Engineering, Northeastern Cniversity, Boston, Mass. 02115
A method for estimating parameter values of linear, unsteady-state distributed models of process behavior in the Laplace domain i s advantageous when closed-form solutions of the equations cannot be obtained in the time domain. In illustrative examples the parameter values of the May model of gas mixing in a fluidized bed are estimated.
E s T I M . k T I x G THE PARAMETERS Of unsteady-state distributed mathematical models so as best to represent actua! data is a problem which occurs frequently in many applications. When the differential form of the model can be solved in closed form in the time domain, the computational work required is usually not excessive. However, closedform solutions often cacnot be obtained and repeated numerical integrations must be performed by computers. Such computations can become prohibitive, since a t least two integrations per parameter arc needed for the determination of a directional derivative a t each step of gradient search, for instance. For linear models with time-invariant parameters it is possible to transform both the model and data into the Laplace domain where the closed-form solution of the model equations can be obtained. The parameters of the model are then estimated with greater ease. The application of the Laplace transform to the problem of estimating model parameters has been used in various cases by others. Chao and Hoelscher (1966) obtained values for the axial Peclet number and mass transfer coefficient in a packed bed from pulse response data. They transformed the model equations into the Laplace domain and obtained analytic solutions for the first three moments of the theoretical response in terms of the model parameters. These parameters were evaluated by numerical computation of the moments from the data. The successful application of this technique depends upon obtaining solutions for the higher moments in terms of the model parameters. This cannot always be accomplished. When such solutions can be obtained, the data in the tailing part of the response curve must be accurate, since relatively small inaccuracies there lead to much larger errors in computing higher moments. Clements (1969) obtained values for the Peclet number of the well-known dispersed plug flow model by performing a least-squared error analysis of the Laplace transformed model and response data along the imaginary axis of the Laplace domain. Parseval’s theorem ensures that the results obtained by this procedure are the same as would be obtained by a least-squared error analysis in the time domain. However, the use of imaginary values of the Laplace variable
Present address, Northeastern University, Boston, Mass.
requires additional mathematical manipulation of the model, since both the real and imaginary parts are involved 111 the fitting procedure. For some models, both the real and imaginary parts of the solution cannot be obtained in practice. Hopkins, Sheppard, and Eisenklam (1969) obtained values of the Peclet number for the dispersed plug flow model using real, positive values of the Laplace variable. Their technique utilizes a logarithmic transformation of the model in the Laplace domain suggested by Ostergaard and Michelsen (1968). This transformation linearizes the model with respect to the parameters and permits them to be evaluated by a linear regression analysis of the data. Their procedure is very effective in this particular case, but in general such linearizing transformations cannot always be found. I n this work it is proposed that the model parameters be estimated by a least-squared error analysis of the transformed response data along the real, positive axis in the Laplace domain. This procedure possesses advantages over those previously discussed. The tailing portion of the response curve is given less weight in the least-squared error analysis than the more accurately determined front section of the response curve. The solution of the model equations in the Laplace domain is used directly without the necessity of obtaining its real and imsgnary parts. Its successful application does not depend upon finding a linearizing transformation of the model solution in the Laplace domain. This method is illustrated by application to the unsteadystate distributed parameter model proposed by hIay to describe gas mixing in a gas-solid fluidized bed. This model was selected because it represents a practical model for which the previously discussed procedures which utilize the Laplace transform for estimating parameter values are not applicable and for which direct estimation in the time domain is very time-consuming. General Procedure
The method proposed here depends upon obtaining a closedform solution of the model equations in the Laplace domain. This solution describes the model response to a known forcing function and is represented in dimensionless form by Equation 1. Ind. Eng. Chem. Fundam., Vol. 9, No. 2, 1970
193
...,
G = G ( S , a i a ~.,
(1)
CY,)
where the ctj are the model parameters. The response data are transformed into the Laplace domain as shown in Equation 2. G d ( S )=
lx’ (-sx)gd(x)dx exp
(2)
Here g d ( X ) represents the dimensionless response curve which is usually available in raw form as a continuous recording. This function, g d ( X ) ,is approximated over segments of the response curve by polynomial expressions and transformed, piecewise, into the Laplace domain. The least-squared error analysis requires that the function, @, be minimized with respect to the parameters of the model, aI,for selected values of S .
analysis arg E , the axial eddy diffusion coefficient, and K , the gas exchange coefficient. The cross-sectional area of the dense-phase compartment, A1, and the volumetric flow rate in the dense compartment, VI, are determined from bed expansion data and from the minimum fluidization velocity, respectively. Model Equations 4 and 5 with their boundary conditions (Equation 6) cannot be solved analytically in closed form. The procedure used to estimate parameter values in the Laplace domain is described below. Model Equations 4 and 5 were transformed with respect to time into the Laplace domain and solved by applying boundary conditions (Equation 6). The solutions in the Laplace domain are straightforward but lengthy; details are given in the appendix. The experimental response curve was closely approximated by linear segments of the form
Y’ 0 The choice of S is not critical, although values of S must be avoided, since the response predicted by any model at S = 0 is insensitive to the selection of the aj.Hopkins, Sheppard, and Eisenklam (1969) recommend values of S between 2 and 5. Their analysis was limited to the dispersed plug flow model and the recommended range of S is specific for that model. In this work it was found that integer values of S from 1 to 10 yielded good results. Any reliable search procedure can be used to minimize the objective function, 9. Levenberg (1944) has discussed a modified gradient search procedure, while Hooke and Jeeves (1958) have described an effective pattern search technique.
afX
+ bi
(7 1
where a i is expressed in terms of the data points ( Y f ’ , X f ) and (Yr+i’,Xi+J by ai = (Yi+I’
- Yi’)/(Xt+l - X f )
(8)
Applying the time-shift theorem, the Laplace transform of Y’ between the limits of X f and Xf+1 can be expressed as
The expression for the transform of the entire RTD data is given by summing over all segments,
P‘ =
Model
=
..N - 1 Pf,,+l‘
=
i=I
May (1959) proposed a model to describe gas mixing in a fluidized bed. The dimensionless model equations, derived from an unsteady-state material balance on a gas tracer, are
-
bY2 bZ
+
(4)
+ (3 (;i:z)2 -
-
-
(Yz- Yl) =
o
(5)
with the boundary conditions Yl(Z,O) = 0 Y,(Z,O) = 0 6atX = 0
(6)
6atX = 0 0 for X > 0 The solution of Equations 4, 5, and 6 yields the response of the fluid bed to a unit impulse input of a gas tracer, Y , where
Y
=
(%)Yl+
(%)YZ
The parameters to be estimated by a least-squared error 194
Ind. Eng. Chem. Fundam., Vol. 9, No. 2, 1970
Further details of this method of numerically transforming the response curve are given by Rooze (1962). The important result to be noted is that the original response curve can be represented in the transformed domain as a function of only the Laplace variable, S . The criterion used to select parameter values is minimization of the sum of the squared deviations between the data values and the model values, both in the Laplace domain. The Laplace variable, S, is a complex number and the fitting procedure may be performed along any path in the plane which does not contain poles of the transformed solution. By Parseval’s theorem, the use of imaginary values of S is equivalent to least squares fitting in the time domain. However, the calculations are simplified considerably if the fitting is performed using real values of S , since the transformed data and model solution will then contain only real values. Paynter (1957) has shown that the poles of systems of the type being considered here lie on the negative, real axis. Hence, in this work, real, positive values of the Laplace variable are used in performing the fit. Integer values of S from 1 to 10 were used for each curve fitting. Although this choice is arbitrary, Jadbabaie (1964) in a related study demonstrated that in the range of S values from 1.0 to 11.0, any extended set of values would give satisfactory results. The use of positive values of S also weights the time domain data so as to give less weight to data points a t larger time values where they are least accurate.
Table I. May Model with Known Parameters Initial Values
3.0 x 10-4 5.0 x 10-5 = 1.0 X lo-* K = 3.0 x 10-3 Fixed parameters V Z / V= 0.28, v2/v
E K E
= =
Final Values
3.528 5.279 3.529 5.279
Actual Parameter Values
3.521 X 5.293 X 3.521 X 5.293 X
X X X X
Sum of Squared Errors
@ = 3.528 X 10-ll @ = 3.538 X lo-”
= 0.53
111ustrat ive Examples The procedure for estimating parameter values in the Laplace domain along the real, positive axis was applied to the May model in three separate cases. First, the model was fitted to data which were generated from the solution of the model equations in the Laplace domain using known parameter values, to determine if the search procedures employed would converge to the correct set of parameter values. Secondly, the model was fitted to response data reported by May (1959). I n this case May obtained a value of parameter E from backmixing data, while the response data were used to estimate a value of parameter K . The purpose was to determine if the same set of parameter values could be obtained by application of the Laplace domain fitting procedure to the response data only. Thirdly, the model was fitted to response data obtained by Williams (1965). The differential form of the May model was integrated in the time domain using the resulting parameter values, to compare the response predicted by the model and the response data in the time domain. Known Parameters. A modified gradient search and a pattern search procedure were tested by fitting the transformed model equations to data generated from selected values of E and K . I n each case the search was initiated from starting parameter values which differed from the true parameter values by a t least an order of magnitude. The results of applying the fitting procedures in the Laplace domain are shown in Table I. The parameter values determined by the search procedures differed from the true values in the third significant figure. The estimated parameter values were found to be unique, although the search was initiated from starting points which differed by two orders of magnitude. Both search procedures yielded essentially identical results. These results show that the search procedures used are capable of minimizing the function, @.
May’s Data. The Laplace domain parameter estimation procedure was applied to gas-phase response data reported by May. May obtained his values for E and K from independent experimental test results. The value of E was obtained from backmixing tests conducted using radioactive solid particles. It was assumed that the eddy diffusion coefficient which characterized the movement of the solid particles also characterized the movement of the gas within the dense phase of the bed. The value of K was obtained from the fluid bed response to a step change of a gas tracer. Only the tailing section of the response curve was utilized in determining K. This was accomplished by numerically solving the model equations for several assumed values of K and comparing the slope of the tailing section of the response curve, plotted on semilog coordinates, to the slopes predicted by the model. The value of K which yielded an agreement between the model and experimental results was selected. May reported a single response curve (May, 1959, Figure 11) taken on a 5-foot diameter fluid bed. In the present work this curve was differentiated and replotted to obtain the bed response to a unit impulse input, so that the model solution developed in the appendix could be applied. This was accomplished by plotting the negative slopes a t selected points along the curve. The curve was truncated a t a value of X T = 2.6 and transformed numerically into the Laplace domain where values for E and K were estimated. The results obtained are shown in Table I1 and compared there to those reported by hfay. The two results are in good agreement, considering the greatly different techniques used in obtaining these parameter estimates. The value of E estimated by fitting in the Laplace domain was somewhat larger than reported by hfay. May reported a “minimumfJ value for E-Le., a value sufficiently large to ensure complete solid mixing in the observed time as measured near the bottom of the bed when the radioactive solids were
Table II. Comparison of Parameters
Laplace domain estimation
Final Values
Initial Values
Method
E K E K
= = = =
2.0 5.0 0.02 0.05
May’s method
E K E
= 6.888
K
= = =
E K
= =
1.216 6.892 1.216 5.0. 0.775b
Sum of Squared Errors
9 = 1.091
x
10-4
@ = 1.091
x
10-4
@ = 4.603 X
Fixed parameters V2/V = 0.424,~v2/v = 0.985 From Figure 10 of May. * From Figure 14 of May. c Based on assumption of 40y0 void fraction in dense phase. a
Ind. Eng. Chem. Fundam., Vol. 9, No. 2, 1970
195
injected a t the top. It is clear from Figure 9 of May’s work that the time required to achieve uniform solid tracer concentration as predicted using May’s value of E is greater than actually observed a t intermediate bed levels. The value of E estimated in the Laplace domain is larger than the “minimum” E and is believed to be a reasonable estimate. The agreement between the values of K is especially good, since the May method of estimating K is very sensitive to experimental error in the tailing section of the response curve. For example, an instrumental or recording error of 2y0 of full scale would shift the slope of the step response curve sufficiently to bring the two estimated values of K into exact agreement. Likewise, an error of 2y0 in plotting May’s data on similog paper to determine the slope could also account for the difference between the values of K obtained. May’s parameters yielded a sum of squared errors in the Laplace domain an order of magnitude larger than obtained using the parameters estimated by the Laplace domain fitting procedure. This suggests that the use of the parameters estimated in the Laplace domain would result in a better comparison of model and data in the time domain for the front sectionof the response curve than the use of May’s parameters. Pulse Test Data. May’s model was used to estimate parameters E and K by fitting it to gas phase response data taken on a 3-inch-diameter fluid bed using radioactive krypton-85 as the tracer. Details of these measurements are discussed by Killiams (1965). The fluidization conditions under which the data were taken are summarized in Table 111. The estimated parameter values are E = 0.0035 and K = 0.00053. The values of E and K are several orders of magnitude smaller than the corresponding values obtained in fitting Xfay’s data. This result is expected for two reasons. First, May showed that E decreases with decreasing bed diameter. Secondly, the 3-inch bed was fluidized a t a much smaller gas velocity than May’s fluid bed. This can be appreciated by referring the ratio of the gas velocity to that required for minimum fluidization. I n the case of the May data, this ratio was 90; for the Williams data, it was only 2.7. The value of K tends to 0 as this ratio tends to 1.0. The estimated parameter values were used in the differential form of the time domain model, Equations 4,5, and 6, which were numerically integrated. The integrated response curve is shown in Figure I, where it is compared to the data. The two results are in good agreement. Although this result does not conclusively establish that the estimated parameters are physically significant and unique, it is a necessary condition that this be the case. It does demonstrate that parameter estimates obtained from fitting in the Laplace domain can lead to good agreement between model and data in the time domain.
Table 111. Fluidization Conditions Solid particles
Microspheroidal cracking catalyst
hlean diameter of particles, microns Gas Pressure, p.s.i.g. Temperature, O F . Gas flow rate, ft.asee.-’ Minimum gas flow rate, ft.3 set.-'
30.7 Air 15 83 7.99 x 10-4 3.0 x 10-4
Numerical errors arising from the transformation of the front section of the response curve into the Laplace domain can be controlled by using as many segments to approximate the curve as desired. Experimental errors in the tailing portion of the response curve are not significant in affecting the values of the estimated parameters. The effect of experimental error in the tailing section of the response curve can be better appreciated b,y reference to a simple example. The dimensionless response curve for a perfectly mixed flow system to a unit impulse input is repre1 sented byG(S) = in the Laplace domain. If the response ~
S+1
curve were truncated a t a dimensionless time value of X r , the transformed response curve would be represented by the
cation is equivalent to assuming thltt the tailing section of the response curve cannot be determined experimentally. The squared error contributed by truncation is given by [G(S) GT(S,XT)I2=
{
[-(’
S+1
1)X’1}2. The contribution to the
sum of squared errors would be of the order of magnitude 10-6 when X T = 3 and 10-8 when X T = 4. Errors of this magnitude would have a negligible effect on the estimated parameter values of the X a y model. In any particular case the error contributed by the tailing section of the response curve can be estimated by truncating the curve a t different points and comparing the results obtained by numerical transformation. The method of fitting in the Laplace domain can also be used in combination with other parameter estimation techniques such as that used by May. In that case the backmixing results could be used to establish a value of E , while the step response curve could be fitted in the Laplace domain to
Discussion
The general approach to parameter estimation through transformation of coordinates so as to estimate model parameter values with greater ease is well known. Logarithmic transformations which linearize the model with respect to the parameters so that linear regression analysis may be used are common in chemical engineering applications. Such transformations are made a t the expense of distorting model and data with the objective of facilitating the estimation procedure. The method of coordinate transformation proposed in this work distorts the data so as to give more weight to the front section of the response curve and less weight to the tailing section where the experimental data are least accurate. 196 Ind. Eng. Chem. Fundam., Vol. 9, No. 2, 1970
Y
5
X Figure 1. model
Comparison of data with RTD computed from the
obtain K . This approach would be useful in the case of multiparameter models where reliable independent data could be used to establish values for individual parameters.
Pz(0,S) = 1 The solutions of Equations 11, 12, and 13 for 2
PI =
Acknowledgment
The authors thank the Systems Research Center, Case Western Reserve university, for financial support in the form of a Graduate Assistantship (to John A. Williams), and the Computer Centers, Northeastern University and Case Western Reserve University, for making available the digital computing facilities used in performing the computations reported here.
where Pi,i
Pa
kl exp (01)
=
+ k~ exp + k3 exp (02)
=
1 are
(P3)
(14)
1, 2, 3, are the roots of the cubic equation,
+ P2(HUS +
J2
- F ) - P [ ( H i s + JI)
+ G ( H 8 + Jd1 (HiHUS2 + JiHUS + JzHiS) 0 =
Nomenclature
A = cross-sectional area, ft2.
c = concentration of tracer, lb,.
ft-3. E = eddy diffusion coefficient, ft.2 sec-I. K = cross-mixing Coefficient, ft.a sec-I. ft-’. L = length of bed, ft. m = number of model parameters N = number of segments used to approximate response .. curve P = number of data points used in least-squared error analysis Q = quantity of tracer input, lb,. s = Laplace variable, dimensionless t = time, sec. v = volumetric flow rate, ft.3sec-I. v = volume, ft3. x = axial distance from bed entrance, f t . X = dimensionless time, vt/V X r = dimensionless truncation point Y = dimensionless concentration, C V / Q 2 = dimensionless axial distance, x / L (Y) = j t h model parameter 6 = impulse function, dimensionless @ = sum of squared errors in Laplace domain
and F
=
VI L EA1 ~
vL2 H -- VE
. -
SUBSCRIPTS 1 = dense-phase compartment of model 2 = void compartment of model i = index of data points = index of model parameters j k = index of Laplace variable
SUPERSCRIPT
’
=
experimental data point
Appendix
hfodel Equations 4 and 5 are transformed with respect to X to yield Equations 11 and 12, respectively.
literature Cited
dPr dr
ilzvL KL +( z SP2 ) + (--)
( 7 2
-
71) =
The transformed boundary conditions are
P,(Z,O)= 0 P,(Z,O)
=
P,(O,S)-
0
y$) 2 j
z=o
= 1
0
(12)
Chao, R., Hoelscher, H E., A.I.Ch.E. J. 12, 271 (1966). Clements, W. C., Jr., Chem. Eng. Sci. 24, 957 (1969). Hooke, R., Jeeves, T. A,, Oper. Res. 6, 881 (1958). Hopkins, M. J., Sheppard, A. J., Eisenklam, P., Chem. Eng. Sci. 2 4 , 1131 (1969). Jadbabaie, hf. J., Ph.D. thesis, Case Institute of Technology, 1964. Levenberg, K., Quart. Appl. Math. 2, 164 (1944). May, W. G., Chem. Eng. Progr. 55, No. 12, 49 (1959). Paynter, H. XI., “On an Analogy between Stochastic Processes and Monotone Dynamic System,” p. 243, “Regelungstechnik: Moderne Theorien und ihre Verwendbarkeit,” Verlag R. Oldenbourg, Munich, Germany, 1957. Ostergaard, K., Michelsen, 31. L., Tripartite Conference, 1968. Rooze, J. W., M. S. thesis, Case Institute of Technology, 1962. Williams, J. A., Ph.D. thesis, Case Institute of Technology, 1965. RECEIVED for review March, 1968 ACCEPTEDFebruary 4, 1970 Ind. Eng. Chem. Fundam., Vol. 9, No. 2, 1970
197
