2175
I n d . E n g . Chem. Res. 1988, 27, 2175-2179
Nonlinear Transformations for Parameter Estimation David M. Espie and Sandro Macchietto* Department of Chemical Engineering and Chemical Technology, Imperial College, London SW7 PBY, United Kingdom
Nonlinear parameter transformations are considered for least-squares estimation in univariate models. A method for the comparison of different transformations is evaluated. When no straightforward transformations are available a priori, a simple family of power transformations is used. Two measures of nonlinearity are used in combination to identify the most suitable transformation, and working constants in the transformations are adjusted to minimize the parameter effects curvature. Several examples are given to illustrate the advantages of this approach. Effective transformations are shown to reduce the parameter effects curvature, the number of iterations required, and the dependency upon good initial guesses of the parameters. 1. Introduction In this paper we consider the use of nonlinear parameter transformations for parameter estimation in univariate models. Ad hoc nonlinear parameter transformations have often been used by experienced users to improve the computational efficiency when fitting models to experimental data (Biegler et al., 1986). Researchers have been hampered, however, by conflicting results over the advantages of using such transformations and by a lack of any methodology to help in choosing a suitable transformation. For example, Agarwal and Brisk (1985) described computational results showing the advantages of reparameterization upon experimental design calculations, in which the experimental points were chosen to minimize the volume of the parameter confidence ellipsoid. Rimensberger and Rippin (1986), however, criticized the results and commented that there should be no effect from the use of nonlinear transformations, since the optimal experiment is independent of the parameter transformation. The arguments of both Agarwal and Brisk and Rimensberger and Rippin also apply in the case of parameter estimation since the position of the global minimum of the sum of squares surface is independent of the parameter transformation. Our discussion of nonlinear transformations will show that both points of view are correct. It will focus on how suitable transformations can be chosen and clarify how they affect the efficiency and robustness of least-squares minimization algorithms whilst confirming the uniqueness of the optima and eliminating doubts about convergence to local minima. Most modern algorithms for solving the nonlinear least-squares problem are based upon linearizing the model about fixed parameter values. The linear approximation assumes that the model may be approximated by its tangent plane and that a uniform coordinate assumption applies. Nonlinear transformations of the parameters affect this uniform coordinate assumption and it has been suggested that a parameterization which reduces the nonlinearity should give improved convergence for existing algorithms and reduce the dependency upon good initial guesses (Ratkowsky, 1983). We therefore must consider measures which quantify the nonlinearity. Using these measures of nonlinearity, we will show how it is possible to choose a rather general transformation of the parameters that can be adapted to each problem. In section 2 there is a brief description of how nonlinearity is quantified, and in section 3 we define a family of transformations that can be used to reduce a measure of
* Author t o whom correspondence
should be addressed.
088S-5885/88/2627-2175$01.50/0
nonlinearity. It the last section, examples are given to compare the performance of transformed and untransformed models using the same nonlinear least-squares algorithm. 2. Measures of Nonlinearity We are concerned with models of the form Y t = f ( x t , @+ )
et
(1)
where yt is a single variable measured at times t = 1, ..., n,; n, is the number of measurement points; xt is a m-dimensional vector of experimental conditions; et is an additive random vector representing measurement noise; and 0 is a p-dimensional vector of time invariant parameters. If the measurement errors have zero mean, then for each experimental setting the conditional expected response is
%(e)= E[Ytl0l = f , ( X t , 0 )
(2)
The sample space is an n-dimensional vector space where a point is defined by the vector y = (yl, y 2 , ..., y,,). The solution locus is a p-dimensional surface consisting of all points in the sample space with coordinates expressible as (3) d e ) = (%(e), ..., ?,,(e)) The linear approximation assumes that (3) can be approximated by its tangent plane about the point 0, and that a uniform coordinate system applies upon q(e). Several measures of nonlinearity have been defined to quantify how closely the tangent plane and uniform coordinate assumptions approximate the solution locus. In particular, Box (1971) used bias in the parameter estimates as a measure of nonlinearity, and Gillis and Ratkowsky (1978) showed that this measure correctly predicted the bias to the right order of magnitude and that the percentage bias provided an indication of which parameters were responsible for the nonlinear effects. Bates and Watts (1980) used the maximum curvature of the solution locus as their measure of nonlinearity. They defined the maximum intrinsic curvature as
rN= maxh (ijhN/fih2)
(4)
and the maximum parameter effects curvature as
rPE= maxh (ijhPE/fih2)
(5)
where h is a unit vector in the parameter space; f i h is (dqt/de)e=e,h;ijh is hT(d2q,/de2)e=eh ; iihN is the component of ijh normal to the tangent plane; and qhPE is the component of q h parallel to the tangent plane. The intrinsic curvature is a property of the solution locus only. In contrast, the parameter effects curvature depends 0 1988 American Chemical Society
2176 Ind. Eng. Chem. Res., Vol. 27, No. 11, 1988
upon the particular parameterization, and rPE can be reduced by a suitable transformation of the parameters. To determine if the nonlinearity is significant a t the a X 100% confidence level, the measures of Bates and Watts (1980) are compared against the distribution
where F@, n,-p, a ) is a F distribution with p and n,-p degrees of freedom. If is less than d(p, n-p, a ) , then the solution locus may be considered to be sufficiently linear at 0,. Similarly, if rPE is less than d@, n-p, a ) ,then the projected parameter lines of 8 are sufficiently parallel and uniformly spaced. There are no corresponding exact limits for the bias measures of Box (1981), but a useful rule of thumb is that if the percentage bias exceeds 1% then this indicates significant nonlinearity. The advantage of the measures of Bates and Watts (1980) is that there is a simple procedure to recalculate the measures of nonlinearity for various transformations (Bates and Watts, 1981). Thus, when searching for a transformation, a large number of these may be tested efficiently. Box’s bias measures can also be easily calculated from the curvature measures. 3. Choosing a Suitable Transformation In specific situations, the form of the transformation may be self-evident or learned from experience. For example, logarithmic transformations of Arrhenius-type kinetic rate constants are usual, while for the HougenWatson-type rate equations a possible general transformation was suggested by Ratkowsky (1985). In this case, the measures of nonlinearity can be readily calculated and should indicate whether a proposed transformation will bring any computational advantages or not. In the majority of cases, however, user supplied transformations will not be available, and instead we propose to use a family of transformations that can be easily adapted to each problem. Nonlinear regression problems often possess “bananashaped” contours. If we express the shape of a valley on this surface as a hyperbola xlx2 = constant
Y1 = X I Y2 = X l X 2
=
where rCzb is the rate of formation of ethane; xczb and xH are mole fractions of ethylene and hydrogen; and all and cyz are constants. The set of unknown parameters is defined as 0 = (Ao, E , al,q).Agarwal and Brisk (1985) used a logarithmic temperature transformation, where
P3
(7)
is useful. The extension of this form that we propose to use is a general power transformation. For example, with three parameters, p2
4. Numerical Experience Several examples are used to show whether transformations which reduce the nonlinearity measures can reduce the computational effort required and can improve the robustness of a typical least-squares algorithm. The computational effort is measured in terms of the number of iterations and function evaluations. Firstly, we will consider the effects of a common transformation and of a transformation that can be easily defined from the form of the nonlinear equation. Secondly, the power transformations will be tested to see if they are useful in general, and they will be compared to some ad hoc transformations for special cases. Finally, an example will be presented where it is not obvious from the form of the equations alone which transformation would be useful. The leastsquares algorithm used in all cases is NL2SOL (Dennis et al., 1981). 4.1. Example from Agarwal and Brisk (1985). This example shows the effect upon the measures of nonlinearity and the computational effort required to obtain a solution, of a commonly used transformation, which from experience is known to be useful for parameter estimation. Results were given in the reference for a sequential experimental design for precise parameter estimation for the hydrogenation of ethylene over a support nickel catalyst. The original rate equation was
p2 = log
then the transformation
p1 = 4
alone or in combination, can be used when a user supplied transformation is not obvious from the form of the model.
OIOZO
p3 = 0,02~0,b
(8)
The flexibility of (8) arises from the fact that the labels 1-3 on the parameters are arbitrary so that any permutation can be used, whilst adjusting the indexes a and b. Here we follow a two-step procedure for selecting the transformation. The most suitable parameters to be included in the transformations are identified by using Box’s bias measure. This measure is calculated for all the parameters, and the three with the highest percentage bias are then included in the transformation. In many cases, logarithmic transformations were found to be invaluable, so these are also included in the list of possible transformations and are compared to the optimal power transformation. Both of these transformations,
=
E R
“I
0 4 = a2 Tf=--T T*
T*=338K
Using the data from a set of 15 designed experiments and a starting point eo= ( 4 w , 10000,0.5, L2), NLZSOL converged in both cases to 9 = (4030.0, 11473, 0.3405, 0.9911). When a transformation is applied, the initial values of the transformed parameters are obtained from 8,. The least-squares solution is then found in the transformed parameter space, and the final values are converted back to the original space by applying the inverse transformations. Table I shows the relative curvatures and the computational effort required for both the original and transformed parameterizations. Because rPE is much larger than d(4, 11,0.95), there is significant parameter effects nonlinearity both before and after the transformation. However, the transformation reduces rPE by 27%, and the number of iterations and function eval-
Ind. Eng. Chem. Res., Vol. 27, No. 11, 1988 2177 Table I. Results for Example 4.1: 6 = (4030, 11 473, 0.3405, 0.9911). I” = 1.4201. d ( 4 , 11, 0.95) = 0.273
iterations function evaluations
22 29
6 7
Table 11. Comparison between Performance of NL2SOL with Different Parameterization8 (Examde 4.2) log el, e2, 82e3-1.6 01, 029 83 0.2067 X los 0.293 X lo6 rPE iterations 251 33 327 53 function evaluations ~~
uations becomes only one-third of those required for the original formulation (Table I). 4.2. Example 8 of Meyer and Roth (1972). In this example, the form of the transformation is obtained by inspection from the original equation. The model used was
The original starting point was eo = (0.02, 4000, 250). By use of the original parameterization, NL2SOL converged to the least-squares solution 0 = (0.005 61, 6181.4, 345.2) in 251 iterations and 327 function evaluations. This is similar to the results quoted by Dennis et al. (1981) (206 iterations, 335 function evaluations) and Al-Baali and Fletcher (1985) (206 iterations and 342 function evaluations), both using NL2SOL. These slight differences are probably due to different machine implementations, for example, different stopping tolerances, which will affect the performance of the code. The curvature measures of the original parameters about the starting point are = 2.5519 and rPE = 0.2067 X lo9. For virtually all the examples that have been published, the intrinsic curvature is much less than the parameter effects curvature (Ratkowsky, 1983). Similarly, for this example, the intrinsic curvature is negligible compared to the parameter effects curvature. Box’s (1971) bias measures for %1, %2, and e3are 3.3 X lo6, 187.7, and 16.1, respectively, and thus the main contribution to the nonlinearity arises from %1 and suggests that a transformation involving would be useful. Taking the logarithms of both sides of (ll),then the following transformation is found by inspection: = log 81
62
=
82
p3
=
82/%31’5
(12)
This reduces rPE to only 0.293 X lo6. Table I1 shows the comparative performance of NL2SOL on this problem solving for both 0 and 8. In this example, the reduction in the computational effort is very significant. 4.3. Example 1 of Meyer and Roth (1972). In the next two examples, the family of power and logarithm transformations suggested in section 3 are compared to transformations which have been suggested in the literature of reducing the measures of nonlinearity. Example 1 of Meyer and Roth (1972) was
The original starting point was eo (10.39,48.83,0.74), and the least-squares solution is 8 = (3.1315, 15.159, = 0.0375, rPE = 12.809, and 0.78006). At this point, d ( 3 , 13, 0.95) = 0.13. Therefore, the intrinsic curvature is negligible, but the parameter effects curvature is unacceptable.
Table 111. Results for Example 4.3: 6 = (3.1315, 15.159, 0.7soi), r N = 0.035, d(3,13,0.95) = 0.13 1feie3, lie3,
e,, e2, e3 e,, e2, e,e3 rPE 12.809 0.264 iterations 13 7 9 function evaluations 19
e2iw3 0.089 5 4
e 2 e ~ 9e2, , ez0,~.90so.9
0.349 7 10
Table IV. Results for Dehydration of Ethanol (Example 4.4)
rPE
12 229 16 30
iterations function evaluations
10.02 9 17
665.3 14 26
The transformation suggested by Bates and Watts (1981) was 01 = %1 P 2 = %2 P 3 = %1%3 (14) This gave a value of rPE= 0.264, which indicates the transformed function is close to being linear at the 95% confidence level. Ratkoswky (1983) suggested
with rPE = 0.089, which is very close to linear. For this example, we calculated a power transformation using Box’s (1971) bias measures to pinpoint the most nonlinear parameters. The percentage biases about 0, were (271, 0.28, 16.86). Thus, most of the nonlinearity is associated with el and e3, and these two parameters were included in the power transformation. The optimum transformation was found with values of the exponents a = 0.9 and b = 0.9 yielding rPE = 0.349. Table I11 shows the comparative computational results using NL2SOL with the initial parameterization and with three different transformations. The transformation of Ratkowsky is the most efficient for reducing both the parameter effects curvature and the number of iterations required. The power transformation and the transformation of Bates and Watts have similar parameter effects curvatures, and both require a similar amount of computation to find the solution. All the transformations converge in fewer iterations than the original problem. 4.4. Dehydration of Ethanol. Kabel and Johanson (1962) gathered experimental results for the vapor-phase dehydration of ethanol. The rate equation proposed was PA2
r=
el
-
pEpW
Keq
(1
+ %,PA + % 3 P w ) 2
(16)
where P A is the partial pressure of diethyl ether, Pw is the partial pressure of water, and PE is the partial pressure of ethanol. Kq was given as 25.2. The initial point used was eo = (0.01, 25.0, 37.5). About this point, = 6.147 and rPE= 1 2 229. The transformation of Ratkowsky (1985) was 01 = 1/%l1I2 @2 = ~ 9 2 / % 1 ~ / ~ 0 3 = %3/%11/2 (17) This reduced rPE to 10.02. Box’s bias measures for %1, and %3 were 12.213, 22.51, and -22.0, respectively; therefore, %2 is the most nonlinear of the parameters and along with e3 was included in a power transformation. The optimum transformation, however, was found to be a
%2,
2178 Ind. Eng. Chem. Res., Vol. 27, No. 11, 1988
that linear scaling of the parameters had no effect upon the zone of nonconvergence.
10.0,
7.6
i
5. Conclusions
0.0
0.00
0.02
0.04
0.06
0.08
0.10
0.12
61 Figure 1. Zone of nonconvergence for example 4.5 (0, = 0.015).
logarithm transformation of 81, which reduced rPE to 665.3. Table IV compares the performance of NLXSOL upon the different transformations. Again transformations that reduced the parameter effects curvature also produced significant reductions in the required number of iterations and function evaluations. 4.5. Ordinary Differential Equation Model. In this example, we will examine a case where ad hoc transformations are not immediately evident. Consider the following reactor model consisting of ordinary differential equations:
91 = 81~1’ yi(0) = 1.0 YAO) = 0.0 9 2 = 81~1- 8 2 . ~ 2+ 6 0 3 3’3 = 8 ~ -28g3 y3(0) = 0.0 (18) It is assumed that only yz is measured; hence, the model
still conforms to the univariate assumption. Since the model function for y is not given explicitly, the partial derivatives of y with respect to the parameters were found by finite differences. A set of nine data points were generated by integrating (18) from t = 0.0 to 10.0, using the solution 8 = (2.0,1.3,0.3) and adding randomly distributed noise with a standard deviation of 0.2. Figure 1 shows a region of the parameter space from which NL2SOL does not converge to the solution. The nonlinearity measures of the power transformations and logarithm transformations were calculated to find the transformation that would offer the largest reduction in nonlinearity. Starting from the point eo= (0.08,6.0,0.01), the percentage biases for all the parameters were similar, so the order of the parameters in the transformation was unimportant and was kept a t 1, 2, 3. A power transformation was found to have the lowest value of rPE, the grid search found the indexes a = -0.9 and b = 0.9, and NL2SOL then found the correct solution for the transformed parameters in eight iterations. To see why a transformation has been successful, when the model with the original parameterization could not be solved, it is necessary to consider the gradients of the objective function with respect to the parameters. NL2SOL failed to converge because the gradient with respect to 0, is much smaller than those with respect to 81 or 83. A t 8, the gradients with respect to B,, 82, and e3 are 10.8, 1.0 X and 0.5, respectively. The response surface is, therefore, poorly conditioned. By use of the optimum transformation, the gradients with respect to pl, &, and p3 are closer in magnitude (10.8, 3.0, 6.5). The nonlinear transformation has altered the sum of squares surface to allow convergence. We also observed
We have described a systematic methodology for utilizing nonlinear transformations in parameter estimation. A power transformation of the parameters was proposed for the case when no obvious transformation was available. This transformation is tuned automatically for individual problems by optimizing the exponents so that the parameter effects curvature rPEis minimized. The parameters to be used in the transformation are identified by comparing the bias of the parameter estimates. Under these conditions a power transformation is effective for either reducing the number of iteractions required or reducing the dependency upon good initial guesses. The results presented show the validity of the argument of both Agarwal and Brisk (1985) and Rimensberger and Rippin (1986). Provided that the nonlinear transformations are properly chosen, then the numerical performance of existing nonlinear least-squares algorithms can be improved. A reduction in the computational effort is consistently obtained when the Bates and Watts parameter effects measure is reduced; however, at the same time the identical final values of the parameters are obtained. Therefore, while the theoretical arguments hold, the numerical experience indicates that there is a definite advantage in using a properly chosen nonlinear transformation. This is quite analogous to scaling or pivoting in the solution of nonlinear equation systems. Nonlinear transformations can also have other less easily quantifiable advantages, such as improving the robustness of the algorithms by enlarging the region of convergence (as in example 4.5), reducing the dependency upon a particularly robust implementation of a parameter estimation algorithm, or increasing the validity of statistical tests based upon the assumption of a linear model. In these cases, even if there is no appreciable improvement in the performance of a particular algorithm, the transformations still serve a useful purpose for the chemical engineering modeler. Problems can arise with the use of nonlinearity measures because the measures of Bates and Watts may be scale dependent. The absolute value of the curvature measure itself is not a consistent indicator of the likely size of the reduction in computational effort achievable when comparing the nonlinearity of different problems; however, the measure is useful for the relative assessment of nonlinear transformations for the same problem. In addition, the measures of Bates and Watts only apply to univariate models, while a larger number of chemical engineering models utilize multivariate data. There is no easy method of extending the use of the univariate measures to multivariate systems and the curvature measures for the multivariate case need to be derived directly. The parameter effects curvature does not apply just to the model but to the model/data set, and joint transformations of the response and parameter space may also be useful for further reducing the parameter effects curvature. Literature Cited Al-Baali, M.; Fletcher, R. “Variational Methods for Nonlinear Least Squares”. J . Oper. Res. SOC. Am. 1985, 36, 405-421. Agarwal, A. K.; Brisk, M. L. “Sequential Experimental Design for Precise Parameter Estimation. 1. Use of Reparameterization”. Ind. Eng. Chem. Process Des. Dev. 1985,24, 203-207. Bates, D. M.; Watts, D. G. “Relative Curvature Measures of Nonlinearity”. J . R. S. S. E. 1980, 42, 1-25. Bates, D. M.; Watts, D. G. ‘Parameter Transformations for Improved Approximate Confidence Regions in Nonlinear Least
I n d . Eng. C h e m . Res. 1988,27, 2179-2182 Squares”. Ann. Stat. 1981,9, 1152-67. Biegler, L. T.; Damiano, J. J.; Blau, G. E. “Nonlinear Parameter Estimation, a Case Study Comparison“. AIChE J. 1986, 32, 29-45. Box, M. J. “Bias in Nonlinear Estimation”. J. R. S. S . B. 1971,33, 171-201. Dennis, J. E.; Gay, D. M.; Welsch, R. E. “An Adaptive Nonlinear Least Squares Algorithm”. ACM Trans. Math. Soft. 1981, 7 , 348-368. Gillis, P. R.; Ratkowsky, D. A. “On the Behaviour of Estimators of the Parameters of Various Yield Density Relationships”. Biometrics 1978,34, 191. Kabel, R. L.; Johanson, L. N. “Modified Damped Least Squares-An
2179
Algorithm for Nonlinear Estimation”. MChE J. 1962,4621-628. Meyer, P.R.; Roth, P. M. “Modified Damped Least Squares-an Algorithm for Nonlinear Estimation”. J. Znst. Math. Appl. 1972, 9,218-233. Ratkowsky, D. A. Nonlinear Regression Modelling; Marcel Dekker: New York, 1983. Ratkowsky, D. A. “A Statistically Suitable General Formulation for Modelling Catalytic Chemical Reactions”. Chem. Eng. Sei. 1985, 40,1623-28.
Received for review July 22, 1987 Revised manuscript received July 12, 1988 Accepted August 3, 1988
Synthesis of Esters of High Molecular Weight. An Analogue of Jojoba Oil. A Statistical Approach Mercedes Martinez, Emilio T o r r a n o , and J o s e Aracil* Chemical Engineering Department, Complutense University, 28040 Madrid, Spain
The development of an ester, which is the analogue of jojoba oil, which can be used as that product is similar t o sperm whale oil and jojoba oil has been made. A process developed in this laboratory for the production of high molecular weight esters uses the esterification of oleic acid and oleyl alcohol to produce oleyl oleate, which can be used as an analogue of jojoba oil. Two full 23 factorial designs at two levels have been used in the study of the synthesis of the analogue of natural oil. The variables selected for study were reaction temperature, initial catalyst concentration, initial alcohol-acid molar ratio, and fixed reaction time and pressure. The results may be explained in terms of the known chemistry of the system. All conclusions are restricted to the experimental range studied. The amount of ester formed is enhanced by increasing temperature. Empirical models which use first-order polynomial expressions were found to represent the amount of ester adequately, while the detection of the curvature effect indicates that second-order polynomials will be required for a better satisfactory description. The commercial quality of the synthesized product is very similar to jojoba oil. The use of jojoba oil has been increased in the last few years, because of jojoba oil’s close resemblance chemically and physically to sperm whale oil. Jojoba oil is obtained from an evergreen bush of the Buxaceae family, Simmodsia chinensis Sheneider, that grows in semidesert areas and yields a nut that contains 50% oil composed mainly of monoesters of the Cmand Cn alcohols and acids, with two double bonds per chain (docosenyl eicosenoate (Cd2)and eicosenyl eicosenoate (C4&)(Miwa, 1971). Commercial planting of the shrub is developing at a growing pace, and large quantities of oil should be available within 10 years. The number and diversity of uses for which jojoba oil has been proposed are great, such as lubricants for highspeed machinery, numerous pharmaceutical uses, cosmetics, hair oils, etc. (Bhatia and Gulati, 1981). There are two problems in using jojoba oil. The first is that we cannot get a great amount of this oil, and the second is that it is not economical. For this reason, the purpose of this work is to study the synthesis, characterization, and properties of an analogue of jojoba oil. The problem was to get cheaper raw materials and develop a very easy synthesis process. As starting products, we used the most available unsaturated acid and alcohol in Spain, oleyl alcohol (C18)and oleic acid (C18). In the synthesis process, we have used the esterification method using an acid catalyst. We have chosen a catalyst that does not react with the double bonds of the molecules, cobalt chloride (Urteaga, 1985). In this paper, emphasis is placed on the use of factorial design of experiments in order to show that the design and statistical analysis of experiments allow us to obtain a
simple but efficient model for industrial control and to reduce the number and cost of experiments. As an illustration of the methodology, the experimental study of the esterification of oleic acid and oleyl alcohol in a batch reactor was undertaken in order to give the phenomenological relation for the amount of ester obtained as a function of initial temperature, catalyst concentration, and alcohol/acid molar ratio. The variables are the most commonly used for modeling esterification reactions. Experimental Section Equipment. The experiment was carried out in a continuous stirred tank reactor (CSTR) of dimensions 500-cm3semispherical bottom, 7 cm high, and 5 cm diameter. The reactor was equipped with stationary baffles attached along the circumference. A marine-type mixing propeller was employed. The impeller speed was set at 1000 rpm. A temperature recorder and controller and a speed controller were provided. The reactor that was immersed in a constant-temperature bath was capable of maintaining the reaction temperature to within f0.1 OC of that desired for the reaction. Materials Used. The commercial catalyst, cobalt chloride by Merck, of constant activity, was used in our experiment. Oleic acid, purity of better than 99.0%, was supplied by Merck and oleic alcohol, purity of better than 98.070,by Henkel Iberica. Analytical Method. The technique to be used had to be capable of monitoring the reaction components. These products were determined by gas chromatography/mass spectrometry (GC/MS) and quantitatively by
0888-5885/88/262~-2179$01.5Q/Q 0 1988 American Chemical Society
