Parameter Estimation for a Polymerization Reactor ... - ACS Publications

3616

Ind. Eng. Chem. Res. 2004, 43, 3616-3631

Parameter Estimation for a Polymerization Reactor Model with a Composite-Step Trust-Region NLP Algorithm Nikhil Arora and Lorenz T. Biegler* Department of Chemical Engineering, Carnegie Mellon University, Pittsburgh, Pennsylvania 15213

We consider an application of the trust-region SQP algorithm described by Arora and Biegler (Comput. Optim. Appl., in press) to the parameter estimation for a polymerization reactor. In particular, we are interested in the robustness of this algorithm in solving ill-conditioned optimization problems with differential algebraic (DAE) models. For this system, it is possible to measure a variety of states, such as production and concentrations. Some of these states contain more information and lead to more reliable parameter estimates than would be obtained if the states were not measured. To minimize sensor costs, it is thus beneficial to select the smallest combination of state measurements that contains the most information. Because the polymer process model is nonlinear, we decide to test the information content of measurements by analyzing properties of the reduced Hessian and constructing joint-confidence regions. This approach is related to the observability analysis (Albuquerque, J. S.; Biegler, L. T. AIChE J. 1996, 42, 2841. Narasimhan, S.; Jordache, C. Data Reconciliation and Gross Error Detection; Gulf Publishing Company: Houston, TX, 2000. Crowe, C. M. Chem. Eng. Sci. 1989, 44, 2909). However, for nonlinear problems, a reliable and efficient parameter estimation algorithm is essential for the measurement selection procedure. Finally, whereas inputs to the process are often held constant in parameter estimation problems, these inputs often contain noise. It is therefore necessary to reconcile inputs along with states in a simultaneous optimization framework. This problem becomes an errors in variables model (EVM) and can be considerably larger and more difficult to solve. In this study, we also consider the related EVM problem and compare results to the standard parameter estimation problem. 1. Introduction We consider the application of NTREST (nonlinear trust-region estimator),1 a recently developed compositestep trust-region algorithm, for parameter estimation of a polymerization model for on-line estimation. Zacca and Ray5 proposed a model for the production of polypropylene in a loop reactor. They showed that the limiting behavior of the loop reactor could be approximated by a CSTR or a sequence of PFRs. Russo and Young6 presented a simplified version of the above model by approximating the reactor as a CSTR and considering only homopolymerization. They also neglected the effect of different kinds of active sites, which was a unique feature of the work by Zacca and Ray. Bindlish7 analyzed some industrial polymerization processes and performed estimation on a copolymerization process. However, he employed a sequential procedure whose performance depends significantly on the differential algebraic (DAE) solver. Sequential methods also suffer from the possibility of failure of the DAE solver if parameters drift too far. Sirohi and Choi8 and Kiparissides et al.9 have also performed parameter estimation on industrial reactions. Again, the focus has been on either a recursive scheme, such as the extended Kalman filter (EKF), or sequential estimation. Sirohi and Choi also compared the EKF and nonlinear programming (NLP) approaches for parameter estimation. They used the EKF to initialize the NLP problem and * To whom correspondence should be addressed. Tel.: 412268-2232. Fax: 412-268-7139. E-mail: [email protected].

implemented it online. Tatiraju and Soroush10 performed state estimation for a polymerization process in which they compared the EKF and a first-principles model solved by a DAE solver. Ramamurthi et al.11 analyzed the NLP formulation for model predictive control (MPC) of a polymerization process. They used a two-level strategy for estimation and control. On the other hand, simultaneous parameter estimation using dynamic polymerization models is not widely reported. The advantage of a simultaneous formulation is that the model is solved only once. The intermediate steps taken by the NLP solver provide infeasible iterates, and so, models need not be converged far from their optimal values. The large size of a discretized DAE estimation problem and the presence of unstable modes in stiff DAEs warrants the use of special decompositions that exploit the structure of the discretized problem, identify unstable modes, and reduce the effect of such modes. Recently, Cervantes and Biegler12 and Biegler, Cervantes, and Wa¨chter13 developed and implemented an efficient DAE decomposition scheme. This scheme selects the basis matrix automatically or allows the user to select a suitable basis. It also does not need to store the basis or the null-space matrices. One often desires to find the model that describes the process best given the measurement data. Moreover, for processing purposes, it is essential to ascertain the minimum set of measurements that are needed to uniquely estimate the model parameters. The goal of this paper is the parameter estimation of a polymerization reactor with possibly incomplete measurement data. In this case, standard estimation tools often fail

10.1021/ie030237e CCC: $27.50 © 2004 American Chemical Society Published on Web 12/31/2003

Ind. Eng. Chem. Res., Vol. 43, No. 14, 2004 3617

to converge and therefore cannot even determine which parameters can be estimated uniquely. This problem is especially important when only a limited set of states can be measured in the process. To treat this problem, we apply a recently developed algorithm1 that converges even in the absence of complete information. This approach allows us to determine the minimal measurement set for parameter estimation. A more general task not considered in this study is the design of experiments and parameter estimation for the discrimination of competing models. Concepts and background for these tasks can be found in Bard23 and Bates and Watts.24 For dynamic reactor models, Bilardello et al.14 distinguished between different DAE models used to describe a batch reactor. They used the volume of the joint confidence region as the criterion for discriminating among models. This is a common practice employed for the design of experiments. Bilardello et al., however, did not discuss the influence of measurements on parameter estimates and the issue of whether all parameter estimates obtained are unique. In addition, Stewart and co-workers15,16 applied a Bayesian approach to discriminate among rival models of single-response and multiresponse processes. In this approach, they obtained relative posterior probabilities of candidate models for their goodness of fit. They demonstrated this strategy on examples drawn from chemical kinetic models. To determine a minimal measurement set (and also to discriminate among competing candidate models), the key link is a parameter estimation strategy that provides convergence to possibly nonunique solutions. This requires a combination of well-formulated models and a robust nonlinear programming algorithm. Therefore, in this study, our goals are five-fold: (1) to analyze the DAE model and formulate an index-1 system from a high-index DAE, (2) to apply NTREST to the simultaneous parameter estimation of a DAE-constrained model of a reaction occurring in a CSTR, (3) to demonstrate robustness of NTREST on ill-conditioned discretized DAE parameter estimation problems, (4) to distinguish between the information contents of different measurements and to develop a way to obtain a minimum set of informative measurements, and (5) to demonstrate a general estimation EVM estimation problem. In this study, we use the code DAE2NLP17 to discretize the DAE system and develop an implementation with NTREST to estimate parameters of a polymerization reactor model. The next section presents the DAE model for the polymerization reactor and its reformulation to an index-1 system. Section 3 then describes the discretization of this model and the formulation of the parameter estimation problem using DAE2NLP. Section 4 provides a description of NTREST for the solution of the parameter estimation problem. Section 5 provides the setting and data generation of this problem, and section 6 then presents the results for various cases of the standard parameter estimation problem and the EVM problem. Section 7 concludes the paper and describes some directions for future work. 2. Process Description A simplified description of the process is available from Russo and Young,6 and a schematic of the reactor used for the process model is shown in Figure 1. The

Figure 1. Schematic of the reactor for the polymerization model.

set of polymerization reactions is simplified to four major irreversible reactions. Feed to the reactor is considered to be monomer (propylene), activated catalyst, and a transfer agent. The reactions are represented as follows kP

A + M 98 P1

Initiation:

kP

Pi + M 98 Pi+1

Propagation:

kT

Pi + T 98 Di + A

Transfer: Termination:

kD

Pi 98 Di

(1) (2) (3) (4)

Here, A is the activated catalyst, M is the monomer, P1 is a polymer chain of unit length, Pi is a growing polymer chain of length i, T is the transfer agent, Di is a dead polymer chain of length i, kP is the propagation rate constant, kT is the transfer rate constant, and kD is the termination rate constant. Russo and Young applied the quasi-steady-state approximation (QSSA) to the concentration of the growing chains. Following the same approach, we can derive an expression for the concentration of live chains of length 1 as

CP1 )

kPCACM kPCM + kTCT + kD

(5)

We can also relate the growing chains of length i to the growing chains of length (i - 1) by a probability R as

CPi ) RCPi-1

(6)

where R is given by

R)

kPCM kPCM + kTCT + kD

(7)

Now, the QSSA allows the following equations to be derived for the zeroth (λ0) and first (λ1) moments of live chains

λ0 ) λ1 )

CP1 1-R CP1

(1 - R)2

(8) (9)

3618 Ind. Eng. Chem. Res., Vol. 43, No. 14, 2004

We also use the following component mole balances

VC˙ M ) QFCFM - QCM - VkPCACM - VkPCMλ0 (monomer) (10) VC˙ A ) Q

F

CFA

- QCA - VkPCACM + VkTCTλ0 (activated catalyst) (11)

VC˙ T ) QFCFT - QCT - VkTCTλ0

(transfer agent) (12)

VΛ˙ 0 ) -QΛ0 + VkTCTλ0 + VkDλ0 (zeroth dead-chain moment) (13) VΛ˙ 1 ) - QΛ1 + VkTCTλ1 + VkDλ1 (first dead-chain moment) (14) where Cj is the molar concentration of species j ∈ {M, A, T}; Λj is the jth dead-chain moment, j ∈ {0, 1}; V is the reactor volume; QF is the input flow rate of material into the reactor; CFj is the input concentration of species j ∈ {M, A, T} in the feed, and Q is the volumetric outflow from the reactor. We have an ordinary differential equation (ODE) for the overall mass balance

VF˘ ) FFQF - FQ

(15)

[

1 0 0 0 1 0 2 -1/(1 - R) -CP1/(1 - R) 1 3 1/(1 - R)2 2CP1/(1 - R) 0 0 0 0 0 0 0 0

Constant volume:

Number-average molecular weight:

Λ1WM Mn ) Λ0 (17)

λ1WM F

0 1 0 0 0 0 0 -WM/F 0 W MQ

0 0 0 0 0 0 0 0 0 0 0 1 0 0

0 0 0 1 0

0 0 0

0 0 0 0 0 0 0 0 1 WMλ1

]

(21)

Observe that the Jacobian 21 is singular. Row 5 has no nonzero elements, algebraic eq 16 does not have any algebraic variables, and the DAE index turns out to be 2. Hence, if we discretize the DAE system, the Jacobian of this discretized system is ill-conditioned. We use the method of Mattsson and So¨derlind18 to convert the DAE system to index 1 as follows: 1. First, we differentiate eq 16 to obtain

F˘ -

C˙ MWM (FM - FS) ) 0 FM

(22)

2. (a) Then, we define dummy variables: Fˆ :) F˘ and ˆ CM :) C˙ M and (b) rewrite eqs 15 and 22 as

where F is the density of material inside the reactor and FF is the density of the feed. We also have the following algebraic equations

CMWM (FM - FS) (16) F ) FS + FM

0 0 0

FFQF FQ + )0 V V

(23)

ˆW C M M (FM - FS) ) 0 FM

(24)

Fˆ Fˆ -

3. Next, we relate C ˆ ˙ M by M to C

ˆ C M -

QFCFM QCM + + kPCACM + kPCMλ0 ) 0 (25) V V

(18)

4. Finally, we drop eq 15. We now have an index-1 system comprising the following DAEs

(19)

VC˙ M ) QFCFM - QCM - VkPCACM - VkPCMλ0 (26)

where FS is the density of the solids inside the reactor, WM is the molecular weight of the monomer, Mn is the number-average molecular weight, xS is the massfraction of the polymer, and wP is the polymer production rate. Equations 5 and 7-19 describe a set of differential and algebraic equations used to represent the polymerization reactor. We assume that the rate constants are defined by the Arrhenius equation

VC˙ A ) QFCFA - QCA - VkPCACM + VkTCTλ0 (27)

Mass fraction of polymer:

xS )

Polymer production rate:

wP ) λ1WMQ

kj ) kj0e-Ej/RT

(20)

where j ∈ {P, T, D}. Here, the activation energies are fixed, and only the parameters kj0 are estimated. 2.1. Index Reduction. Before we can use a DAE system for estimation, we have to ensure that it is of index 1. Algebraic equations in high-index systems have singular Jacobians in the algebraic variables and thus lead to discretized DAE systems that are poorly conditioned. Here, we arrange the eight algebraic variables in the order y :) {CP1, R, λ0, λ1, Mn, xS, wP, Q} and calculate the Jacobian of the algebraic equations 5, 7-9, and 16-19

VC˙ T ) QFCFT - QCT - VkTCTλ0

(28)

VΛ˙ 0 ) -QΛ0 + VkTCTλ0 + VkDλ0

(29)

VΛ˙ 1 ) -QΛ1 + VkTCTλ1 + VkDλ1

(30)

CP 1 R-

kPCACM )0 kPCM + kTCT + kD

kPCM )0 kPCM + kTCT + kD λ0 λ1 Fˆ -

CP1 1-R

)0

CP1 (1 - R)2

)0

FFQF FQ + )0 V V

(31) (32)

(33) (34) (35)


Figure 2. Radau collocation with three collocation points/element.

Fˆ ˆ C M -

C ˆW M M (FM - FS) ) 0 FM

(36)

QFCFM QCM + + kPCACM + kPCMλ0 ) 0 (37) V V Mn -

Λ1WM )0 Λ0

(38)

xS -

λ1WM )0 F

(39)

wP - λ1WMQ ) 0 F - FS -

CMWM (FM - FS) ) 0 FM

(40) (41)

ncol

z(t) :) zi + hi

Ωki [(t - ti)/hi]z′ki ∑ k)1

Ω′ki [(tij - ti)/hi] ) δij

(45) (46)

ncol

Observe that eqs 31-34 trivially relate one algebraic variable to some differential variables, and we can calculate CP1, R, λ0, and λ1 and use these calculated values in the remaining equations. 3. Formulation of Parameter Estimation Problem 3.1. Collocation on Finite Elements. The dynamic model of the polymerization reactor is discretized using orthogonal collocation on finite elements. We employ Radau quadrature and divide time into a discrete number of finite intervals with a certain number of intermediate points in each interval. The intervals are referred to as finite elements, and the intermediate points within each finite element as collocation points. Within each interval, we approximate the states and controls (or inputs) by polynomials that depend on the number of collocation points. We refer to Figure 2 to illustrate the procedure. Consider a DAE of the form

F(z′(t),z(t),y(t),u(t),θ) ) 0

parameters, and z′(t) represents the time derivatives of z(t). If z(t) ∈ R|z| and y(t) ∈ R|y|, then F(‚) ∪ G(‚) ∈ R|z|+|y|. We discretize time into a set of finite elements: t0 < t1 < ... < tn, t ∈ [t0, tne]. Within each element i, consider intermediate points tji, j ∈ {1, ..., ncol}, tji ∈ (ti, ti+1], ∀j, ) ti+1, for Radau quadrature.19 To derive the tncol i polynomial approximation of the differential variables, we employ the monomial basis representation.20 Here, we approximate the differential variable at time t ∈ (ti, ti+1] as

Ωki (1) ) 1 ∑ k)1

where we have used the compact notation zi :) z(ti), z′ki :) z′i(tik), and hi :) (ti+1 - ti). Ω′ki is a Lagrange polynomial of degree (ncol - 1), evaluated in the finite element i and at the collocation point k. We also require the differential variables to be continuous across element boundaries, and for this condition, we impose continuity constraints ncol

z(ti+1) ) z(ti) + hi·

(43)

z(0) ) z0

(44)

where z(t) is the set of differential variables, y(t) is the set of algebraic variables, u(t) is the set of inputs or control variables, θ is the set of time-independent

Ωik(1)z′ki ∑ k)1

(48)

at the boundary of each finite element. We approximate algebraic and control variables in a particular element by Lagrange polynomials of degree (ncol - 1) ncol

y(t) )

(42)

G(z(t),y(t),u(t),θ) ) 0

(47)

Ω′ki [(t - ti)/hi]yki ∑ k)1

(49)

ncol

u(t) )

Ω′ki [(t - ti)/hi]uki ∑ k)1

(50)

where we have again used the compact notation for y and u. We allow the algebraic and control variables to be discontinuous across elements, and we define the discretized DAE at each collocation point as


∀i ∈ {1, ..., ne}, j ∈ {1, ..., ncol} F(z′ji,zi,yji,uji,θ) ) 0

(51)

G(z′ji,zi,yji,uji,θ) ) 0

(52)

ncol

zi+1 - zi - hi

Ωki (1)z′ki ) 0 ∑ k)1

z0 ) z0

(53) (54)

which is a set of nonlinear algebraic equations. Depending on the fineness of the discretization and the number of DAE equations, the discretized system can become a large-scale problem. A desire to obtain accurate state and control profiles might lead the user to select fine discretization or an uneven placement of elements. For such cases, exploiting the structure of eqs 51-53 can save computational effort and reduce storage. Wa¨chter17 has developed the DAE2NLP code that exploits the structure of the Jacobian of eqs 51-53 and calculates relevant quantities efficiently. The code essentially uses a discretization such as COLDAE21 and uses the development of Cervantes12 to generate the basis matrix and subsequent computations. The DAE constrained parameter estimation problem is now written as

Figure 3. Confidence region between two parameters.

and y represents the dependent variables. Without loss of generality, we assume that yi ∈ Rny, vi ∈ Rnv, θ ∈ Rnp, and h: Rnv X Rnp f Rny. For simplicity, we assume that there are no errors in the independent variables. Assuming that errors in the measurements of the dependent variables, yi, are normally distributed according to N(0, Vy), we can set up the regression problem to estimate parameters θ as

nm

min Φ(z,y,u,θ) :) z′,z,y,u,θ

T -1 M {[z(tl) - zM ∑ l ] Vz [z(tl) - zl ] + l)1

-1 M T [y(tl) - yM l ] Vy [y(tl) - yl ] +

(55)

G(z′ji,zi,yji,uji,θ) ) 0 Ωki (1)z′ki ) 0 ∑ k)1

(57)

We can use the methods of Bard23 and Bates and Watts24 to obtain the covariance matrix of the parameters, Vθ, as well as the quantity

(58)

distributed according to a χ2 distribution with np degrees of freedom. We now obtain an ellipsoid of the form (Figure 3)

z0 ) z0 where the subscript l denotes the particular times at which measurements are taken. We have developed an interface of NTREST with DAE2NLP to solve dynamic optimization problems that generate the results presented in this section. In addition, we obtain first and second derivatives of the objective function and the constraints using ADOL-C.22 3.2. Parametric Inference. Point estimations of the parameters are obtained from the solution of problem 55, and statistical inferences can be made around these points. This can be accomplished by obtaining confidence intervals and joint confidence regions. Here, we apply an inexpensive method to obtain approximate ellipsoidal confidence regions. Ellipsoidal regions are exact only in the case of linear process models. On the other hand, even with nonlinear models, they provide useful qualitative inference measures for the fitted parameters.24 Consider a model of the form

hi(vi,yi,θ) ) 0,i ∈ {1, ..., nm}

-1 M (yM ∑ i - yi)Vy (yi - yi) i)1

Θ :) (θ - θ*)TVθ-1(θ - θ*)

ncol

zi+1 - zi - hi

θ

s.t. hi(vi,yi,θ) ) 0, i ∈ {1, ..., nm}

T -1 M [u(tl) - uM l ] Vu [u(tl) - ul ]}

s.t. F(z′ji,zi,yji,uji,θ) ) 0

nm

min

(56)

where i is the measurement, v represents the independent variables (similar to process inputs or controls),

(θ - θ*)TVθ-1(θ - θ*) e χnp2(R)

(59)

where R is the confidence level, which we set to 95% in this study. We illustrate the ellipsoid in expression 59 in two dimensions using Figure 3. Here, the principal axes of the ellipsoid are denoted by the lines AB and CD. The principal axes indicate the smallest and largest variations in linear combinations of the parameters that lead to equivalent data fitting. Clearly, the smaller an ellipsoid, the better the parameter estimation. 4. Description of NTREST NTREST is based on a successive quadratic programming (SQP) algorithm designed to solve equalityconstrained parameter estimation problems of the form

min f(θ,y) {θ,y}

s.t. c(θ,y) ) 0 L

(60) U

θ eθeθ

where f(θ,y) is an objective function often derived from


a statistical basis; xT ) [θT, yT] is the total set of variables; θ is the set of parameters; and θL and θU are the lower and upper bounds, respectively, on the parameters. Here, we assume that f: Rn f R, c: Rn f Rm, θ ∈ R(n-m), y ∈ Rm, and n > m. The functions f and c are assumed to be at least twice continuously differentiable. Observe that we have bounds only on the parameters. Problem 60 is solved by successive quadratic programming (SQP) by repeatedly solving the following quadratic program for the kth iteration

where wk ) ZkTWkYkpY, Bk ) ZkTWkZk, d ) YkpY + ZkpZ, pY ∈ Rm, and pZ ∈ R(n-m). To improve convergence properties, the above subproblems are extended to include trust regions for pY and pZ. The resulting composite-step trust-region algorithm consists of the trust-region quasi-normal and trust-region tangential subproblems in NTREST.

Trust-region quasi-normal subproblem min ||ck + AkYkpY||2 pY

1 min ∇xf(xk)Td + dTWkd d 2 s.t. c(xk) + A(xk)d ) 0

(61)

θL e θk + dθ e θU Here, d ∈ Rn, dT ) [dθT, dyT], dθ ∈ R(n-m), and dy ∈ Rm. In what follows, we use the following notation for brevity: fk ) f(xk), gk ) ∇xf(xk) ∈ Rn, ck ) c(xk), Ak ) A(xk) ) ∇xckT ∈ Rm×n, and Wk ∈ Rn×n is the Hessian of the Lagrangian, Lk ) fk + ckT λk ∈ Rn. We define the null space of the Jacobian Ak as the matrix Zk. Here, Zk ∈ Rn×(n-m), and AkZk ) 0. Matrix Yk ∈ Rn×m is calculated so that the matrix [ZkYk] is nonsingular. By suitably arranging the columns of the Jacobian such that

Ak ) [Nk Ck], Ck ∈ Rm×m non-singular, and Nk ∈ R

m×(n-m)

(62)

we can calculate Yk and Zk based on the coordinate decomposition,25 where we have

Zk )

[

I(n-m) -Ck-1Nk

]

[]

0 Yk ) I m

and

(63)

From the first-order necessary conditions of 61, we obtain T

T

Zk Wkd ) -Zk gk

(64)

YkTWkd + YkTAkTλk ) -YkTgk

(65)

Akd ) -ck

(66)

By splitting the step taken by the QP into two parts using the coordinate basis d ) YkpY + ZkpZ, we obtain

[

][ ] [ ]

ZkTWkZk ZkTWkZk 0 pZ ZkTgk YkTWkZk YkTWkYk YkTAkT pY ) - YkTgk λk AkYk ck 0 0

(67)

which leads to the quasi-normal and tangential subproblems as follows

Quasi-normal subproblem ck + AkYkpY ) 0

(68)

Tangential subproblem 1 min (ZkTgk + wk)TpZ + pZT BkpZ pz 2 pLZ e pZ e pU Z

(70)

s.t. ||pY|| e δ1 This subproblem is a large NLP (in Rm) with a quadratic objective and a single quadratic constraint. We assume Ak to be full row rank so that this problem is strictly convex. Because of the large size of this problem, we prefer to solve it inexactly with the aid of sparse matrix factorizations within a dogleg method.

Trust-region tangential subproblem 1 min (ZkTgk + wk)TpZ + pZTBkpZ pZ 2 s.t. || pZ|| e δ2

(71)

θL e θk + pZ e θU This subproblem is a small NLP (in R(n-m)) with a quadratic objective, a single quadratic constraint, and bounds. For this problem, we apply a classical trustregion algorithm extended to handle simple bounds.1 In particular, an important feature of this algorithm lies in computing descent directions when Bk is singular or indefinite. 4.1. Lagrange Multipliers and Second-Order Correction. We can apply a Newton step to calculate the step d and Lagrange multipliers λk of QP 61 at iteration k to obtain

[

][ ] [ ]

YkTWkYk YkTWkZk YkTAkT pY YkTgk T T pZ ) - Z Tg Zk WkYk Zk WkZk 0 k k λk ck AkYk 0 0

If we ignore the quadratic terms YkTWkYk and YkTWkZk in the first equation, we can obtain Lagrange multiplier estimates using

λk ) -(YkTAkT)-1YkTgk

(73)

Here, Lagrange multipliers play a subordinate role in our algorithm and are used only to calculate the exact Hessian of the Lagrangian. NTREST uses the nondifferentiable l2 merit function

φ(x,µ) ) f(x) + µ||c(x)||

(74)

with µ being the penalty parameter. This merit function is susceptible to the Maratos effect, which we counter by calculating a second-order correction step that is added to the quasi-normal step. We employ the cheap calculation

pSY ) -Ck-1c(xk + d) (69)

(72)

(75)

This ensures a move toward feasibility of the constraints.


4.2. Predicted and Actual Reductions, TrustRegion Update, and Penalty Parameter. After calculating the quasi-normal and tangential steps, we assess the quality of the overall step by comparing the reduction, if any, in the merit function to the reduction in the quadratic models. The predicted reduction in the quadratic models contains the predicted reduction in the model for the tangential step (mT)

1 mkT ) -(ZkTgk + wk)TpZ - pZTBkpZ 2

(76)

and the predicted reduction in the model for the quasinormal step (mN)

mN k ) ||ck|| - ||ck + AkYkpY||

(77)

The overall predicted reduction is now

1 T T T predk ) mkT + µ+mN k - gk YkpY - pY Yk WkYkpY 2 (78)

Table 1. Nominal Parameter and Input Values quantity

value

units

kP 0 kT0 kD0 EPA ETA ED A CFM CFA CFT QF FF T FM FS WM V

4.9700 × 4.4000 × 106 7.9200 × 103 1.2000 × 104 1.2000 × 104 1.2000 × 104 8.5714 1.0000 × 10-4 2.7200 × 10-3 1.1554 × 101 3.6001 × 102 3.4300 × 102 3.6000 × 102 8.4000 × 102 4.2000 × 101 5.0830 × 104

L/(gmol s) L/(gmol s) s-1 cal/gmol cal/gmol cal/gmol gmol/L gmol/L gmol/L L/s g/L K g/L g/L g/gmol L

107

where µ+ is the penalty parameter in eq 75 updated to ensure that predk is positive and that µ is positive and monotonically increasing as the solution is approached.1 The actual reduction is just the difference in the merit function at two successive iterates, i.e.

aredk ) f(xk) - f(xk+YkpY+ZkpZ) + µ+(||c(xk)|| - ||c(xk+YkpY+ZkpZ)||) (79) Here, pY might include the second-order correction. The quality of the step is now assessed, and the radius of the trust region altered. Also, we define

Fk )

aredk predk

(80)

as a gauge for increasing or reducing the size of the trust region. Large values of Fk allow the trust region to double, whereas small values of Fk require this region to be cut in half. Details on the adjustment of the trust region are described in Arora and Biegler.1 4.3 Overall Algorithm We are now ready to state the complete algorithm. 0. Set all tolerances and constants, as well as x1, µ1, and δ1. For k ) 1, ... 1. At xk, evaluate fk, ck, gk, Ak, Yk, and Zk. 2. Calculate λk from eq 74, Wk, and Bk ) ZkTWkZk. 3. Calculate YkpY using the dogleg method and set wk ) ZTk WkYpY. 4. Solve for pZ and hence ZkpZ. 5. Update the penalty parameter (µ+) using eq 79 and calculate the merit function fk + µ+||ck||. 6. Calculate ared and pred from eqs 80 and 79, respectively. 7. Check if a second-order correction is needed. Define d ) YkpY + ZkpZ + YkpSY if the secondorder correction is needed or d ) YkpY + ZkpZ if not. Recalculate ared if the second-order correction is needed. 8. Evaluate the step and update the trust-region parameters, δ1 and δ2, according to Fk.1 9. If converged, STOP. Else, if (δ1 < δmin), STOP. 10. Calculate xk+1 ) xk + d. End For

Figure 4. Trended input profiles used for simulation as a function of time (units of seconds on the x axis).

With a description of NTREST in mind, we now briefly summarize how to analyze parameter estimates using joint confidence regions. 5. Simulation and Data Generation In lieu of measurement data for this system, we perform simulations using nominal values of inputs and parameters to generate these measurements. We then add random noise to the simulation results and the inputs and use the noisy data for estimation. We use the nominal parameter and input values (Table 1) given in Russo and Young.6 We perform simulations assuming two types of input profiles. First, we assume constant inputs set to their nominal values over the entire duration of the simulation. Second, we assume trended input profiles for CFM, CFA, and T. Here, we ramp up CFM from t ) 0 for a period of 300 s and then hold it constant. We start with the nominal value of CFA, step it down at t ) 400.00 s, hold it constant, and step it up again at t ) 5940.00 s. Finally, we ramp up T throughout the duration of the simulation. Figure 4 depicts these trended input profiles. Simulations are performed over a period of 5.4 × 104 s. For these simulations, we assume that the reactor

Ind. Eng. Chem. Res., Vol. 43, No. 14, 2004 3623 Table 2. Results of Parameter Estimation When Inputs to the System Are Constant measurement set CM, CT, Mn CM, CT, wP CM, CT, Q CM, CT, Mn, wP CM, CT, wP, xS CM, CT, Mn, xS Mn, Q, xS, wP CM, CT, Mn wP, Q, xS

eigenvalues of reduced Hessian 10-17,

10-9,

10-3

5.18 × 2.90 × 1.10 × 2.59 × 10-13, 5.88 × 10-10, 4.97 × 10-3 5.11 × 10-17, 7.25 × 10-12, 1.44 × 10-3 3.74 × 10-11, 4.38 × 10-9, 5.43 × 10-3 2.50 × 10-13, 8.37 × 10-10, 6.36 × 10-3 8.69 × 10-12, 4.57 × 10-9, 2.23 × 10-3 5.88 × 10-11, 6.60 × 10-9, 7.89 × 10-3 6.14 × 10-11, 6.93 × 10-9, 8.15 × 10-3

initially contains only monomer and transfer agent with CM(0) ) 8.5714 gmol/L and CT(0) ) 2.72 × 10-3 gmol/L. We also assume that CA(0) ) 1.0 × 10-10 gmol/L, Λ0(0) ) 1.0 × 10-6 gmol/L, and Λ1(0) ) 1.0 × 10-6 gmol/ L. To the simulated results of states S :) {CM, CT, Mn, wP, Q, xS}, we add noise sampled from the normal distribution N(0, σj), where σj ) 0.1 × O(jsim), j ∈ S, and the subscript sim refers to the simulated result. We also add noise to the inputs using N(0, 0.005 × |j|), j ∈ {CFA, CFM, CFT, QF, FF, T}. If we also estimate some inputs along with parameters, we add noise from N(0, 0.05 × |j|), j ∈ {CFA, CFM, CFT, QF, FF, T}. For the standard parameter estimation problem, we have to fix the inputs at element boundaries. This is easy to do because we do not add any noise to the input profiles. However, when the input profiles have noise, element placement for parameter estimation is independent of element placement for simulation. Hence, noisy input measurements that are available at certain specific times might not coincide with element boundaries. We thus have to locate two input measurements that are closest to an element boundary, on either side of it, and then interpolate to obtain an input at the element boundary. Inputs at collocation points within an element are obtained by linearly interpolating inputs at the boundaries of that element. 6. Parameter Estimation We present parameter estimation results for two types of data generated by constant and trended inputs. We first consider parameter estimation when all inputs are fixed at element boundaries. In all cases, we use two collocation points. We analyze the quality of estimation with plots of predicted vs measured states, eigenvalues of the reduced Hessian at optimality, and joint confidence regions. We measure up to six states defined by S. The above analysis also identifies the most informative measurements. To prevent the parameter estimates from being estimated too far from their true values, we add the following bounds

3 × 107 e kP0 e 7 × 107 L/(gmol s) 3 × 106 e kT0 e 7 × 106 L/(gmol s) 6 × 103 e kD0 e 7 × 104 L/s 6.1. Constant Inputs. First, we consider parameter estimation when inputs to the process are held constant. Here, we discretize the time interval t ∈ [0, 54 000] into 30 finite elements placed at varying distances from each other. We place the first 10 elements at a separation of 25 s, the next 15 elements at a distance of 150 s, and the last 5 elements at a distance of 10 300 s. We select

parameter estimates {kP, kT, kD} 6.26 × 107, 5.56 × 106, 1 × 104 4.61 × 107, 4.69 × 106, 7.16 × 103 5.36 × 107, 5 × 106, 8.53 × 103 4.99 × 107, 4.43 × 106, 7.83 × 103 4.63 × 107, 4.64 × 106, 7.26 × 103 5.06 × 107, 4.49 × 106, 7.94 × 103 4.96 × 107, 4.37 × 106, 7.90 × 103 4.96 × 107, 4.38 × 106, 7.89 × 103

Figure 5. (CM, CT, wP) reconciled.

such a scheme because the simulations seem to reach steady state or close to steady state fairly quickly and most of the transient information is located within the first 2000 s. We summarize the eigenvalues of the reduced Hessian and parameter estimates obtained by using the different measurement combinations in Table 2. Small eigenvalues of the reduced Hessian indicate that certain combinations of parameters are essentially nonunique at the solution of NLP 55. Small eigenvalues also lead to an almost singular reduced Hessian, indicating that classical NLP algorithms that invert a Hessian would fail or lead to extremely long steps that would have to be significantly cut by line searches. The eigenvalues of the reduced Hessian indicate that some linear combinations of the parameters contain significantly more uncertainty than others. The linear combinations can be derived from the eigenvectors of the reduced Hessian. Because we have used good initial guesses for the parameter values, we obtain parameter estimates that are close to the nominal values reported in Table 1. However, the eigenvalues of the reduced Hessian vary considerably depending on which measurements are reconciled. The eigenvalues increase as we add more measurements to be reconciled and tend to settle down around O(10-11, 10-9, 10-3). The measurement wP seems to be very informative, as the eigenvalues are larger when wP is reconciled than when it is not. For example, the eigenvalues are significantly larger when (CM, CT, Mn, wP) are reconciled than when (CM, CT, Mn) or (CM, CT, Mn, xS) are reconciled. Now consider the quality of the fits presented in Figures 5-8. Also, we obtain qualitatively similar state profiles when we reconcile (CM, CT, Mn), (CM, CT, Mn, wP), (CM, CT, Mn, xS), (Mn, wP, Q, xS), and (CM, CT, Mn,


Figure 6. (CM, CT, Q) reconciled.

Figure 7. (CM, CT, Mn, wP) reconciled.

Figure 8. (CM, CT, wP, xS) reconciled.

wP, Q, xS). Hence, we consider only one figure, Figure 7, as representative of these cases and refer to this figure whenever needed.

Figure 9. Confidence region between kP0 and kT0 when (CM, CT, Mn) reconciled.

First, we compare the estimation with (CM, CT, Mn) to that with (CM, CT, wP) in Figure 5. If we reconcile only (CM, CT, wP), we are unable to fit Mn wellsthe fitted values are underestimated (see Figure 5). The other profiles are fitted similarly. Now, if we compare (CM, CT, wP) and (CM, CT, Mn, wP) from Figures 5 and 7 with Figure 6, where (CM, CT, Q) are reconciled, we observe that neither wP nor xS is fitted well in Figure 6. This case, reconciling (CM, CT, Q), also gives the smallest set of eigenvalues of the reduced Hessian. Comparing these three cases to that shown in Figure 7, we observe that all profiles are fitted well in Figure 7, where four measurements are reconciled. Figure 7 also gives better data fits than does reconciliation of the measurements (CM, CT, Mn, xS) and (CM, CT, wP, xS). Also, the reconciliation of these four algebraic variables (Mn, wP, Q, xS) is better than the reconciliation of any other combination of four measurements that we have tried; the quality of the fit in this case resembles that obtained by reconciling six measurements in Figure 7. This indicates that these four measurements contain sufficient information to give reliable parameter estimates and good fitted values. We now analyze the confidence regions corresponding to the above cases. We first consider the confidence regions between kP0 and kT0. Figures 9 and 10 correspond to the reconciliations of (CM, CT, Mn) and (CM, CT, wP), respectively, and the latter figure shows a greater uncertainty in kP0 and a lower uncertainty in kT0 than the former. The overall uncertainly in kP0 is around 1.6%, and that in kT0 is 1.5% for Figure 9; the corresponding values are 1.2 and 2.0%, respectively, for Figure 10. Figure 11 shows a very large uncertainty in kT0 and a somewhat reduced uncertainty in kP0 for the reconciliation of (CM, CT, Q). However, the increased uncertainty in kT0 indicates that this parameter is essentially nonunique. When we reconcile different combinations of four measurements, as reported in Table 2, we obtain the smallest confidence regions and lowest parametric uncertainty when we reconcile both Mn and wP, as shown in Figure 12. Finally, we compare the confidence regions produced by reconciling four algebraic variables and six measurements (Figures 13 and 14, respectively). The ellipse produced by reconciling six measurements (Figure 14) is only slightly


Figure 10. Confidence region between kP0 and kT0 when (CM, CT, wP) reconciled.

Figure 13. Confidence region between kP0 and kT0 when (Mn, wP, Q, xS) reconciled.

Figure 11. Confidence region between kP0 and kT0 when (CM, CT, Q) reconciled.

Figure 14. Confidence region between kP0 and kT0 when CM, CT, Mn, wP, Q, xS) reconciled.

Figure 12. Confidence region between kP0 and kT0 when (CM, CT, Mn, wP) reconciled.

smaller than the ellipse in Figure 13. This reiterates that only four measurements are required to provide sufficient information for parameter estimation.

The fitted profile for wP in Figures 5-8 is not very good, especially toward the end, suggesting that the poor-quality fitting might be due to coarse discretization. We thus performed parameter estimation with 60 finite elements to test this hypothesis (Figure 15). Here, we reconciled (Mn, wP, Q, xS). The profile for wP became worse, although the profiles for the other states fit their true values very closely. To illustrate how the uncertainty in parameters carries over to combinations of parameters other than kP0 and kT0, we also compared confidence regions between kP0 and kD0 and between kD0 and kT0 for cases 1 and 8 in Table 2. The results are presented in Figures 16 and 17 for the confidence region between kD0 and kT0 and in Figures 18 and 19 for the confidence region between kP0 and kD0. We compare the information content of measurements (CM, CT, Mn) and (CM, CT, Mn, wP, Q, xS). Again, a similar trend is observed in the sizes of the confidence regions. The sizes are significantly reduced when more measurements are reconciled. 6.2. Trended Inputs. We also performed parameter estimations assuming some inputs to the reactor to have trended input profiles. The results of these estimations are presented in Table 3.

3626 Ind. Eng. Chem. Res., Vol. 43, No. 14, 2004 Table 3. Results of Parameter Estimation When Inputs to the System Are Trended measurement set CM, CT, Mn CM, CT, wP CM, CT, Mn, wP CM, CT, wP, xS Mn, Q, xS, wP CM, CT, Mn wP, Q, xS

eigenvalues of reduced Hessian 10-17,

10-9,

10-4

5.18 × 2.90 × 8.58 × 9.11 × 10-13, 1.61 × 10-9, 6.12 × 10-3 5.56 × 10-11, 5.53 × 10-9, 6.48 × 10-3 9.64 × 10-13, 1.77 × 10-9, 7.30 × 10-3 6.46 × 10-11, 6.90 × 10-9, 8.34 × 10-3 6.76 × 10-11, 7.37 × 10-9, 8.78 × 10-3

parameter estimates {kP, kT, kD} 6.26 × 107, 5.56 × 106, 1 × 104 4.61 × 107, 4.56 × 106, 7.59 × 103 4.98 × 107, 4.39 × 106, 8.08 × 103 4.59 × 107, 4.60 × 106, 7.46 × 103 5.01 × 107, 4.44 × 106, 8.01 × 103 5.01 × 107, 4.44 × 106, 8.01 × 103

Figure 15. (Mn, xP, Q, xS) reconciled with 60 finite elements.

Figure 17. Confidence region between kD0 and kT0 when (CM, CT, Mn, wP, Q, xS) are reconciled.

Figure 16. Confidence region between kD0 and kT0 when (CM, CT, Mn) are reconciled.

Figure 18. Confidence region between kP0 and kD0 when (CM, CT, Mn) are reconciled.

Observe that the eigenvalues of the reduced Hessian follow the same trend as in Table 2 for constant inputs. The addition of Mn and wP to the reconciled measurements significantly increases the smallest eigenvalue. The final parameter estimates are also close to the nominal values. We then analyze the fitted values in Figures 20-23. To generate these results, we used 50 finite elements and included a small element with end points closely covering the discontinuity in the measurement profiles. Again, qualitatively similar fitted values are obtained when (CM, CT, Mn, wP), (Mn, wP, Q, xS), or (CM, CT, Mn, wP, Q, xS) were reconciled, and we select Figure 22 to represent all of these cases. We now consider Figures 20 and 21. When wP is not reconciled, as in Figure 20, we observe that the fitted values of wP

and xS are underestimated significantly. When Mn is not reconciled, as in Figure 21, the fitted values of Mn are also underestimated. Again, reconciling four measurements with Mn and wP included as two of the reconciled measurements produces better data fits than reconciling four measurements with either of them not included, as can be observed in Figures 22 and 23. Finally, reconciling four algebraic variables seems to provide sufficient information. We then analyze the confidence regions in Figures 24-27. Consider first Figures 24 and 25. Here, when wP is reconciled, we obtain a smaller confidence region than when it is not. The smaller confidence region in Figure 25 indicates uncertainties of 0.08% in kP0 and 1% is kT0. We compare different combinations of four


Figure 19. Confidence region between kP0 and kD0 when (CM, CT, Mn, wP, Q, xS) are reconciled.

Figure 22. (CM, CT, Mn, wP) reconciled: trended inputs.

Figure 20. (CM, CT, Mn) reconciled: trended inputs.

Figure 23. (CM, CT, wP, xS) reconciled: trended inputs.

Figure 21. (CM, CT, wP) reconciled: trended inputs.

Figure 24. Confidence region between kP0 and kT0 when (CM, CT, Mn) are reconciled: trended inputs.

measurements and observe that reconciling Mn and wP leads to smaller confidence regions. Moreover, reconciling the four algebraic variables produces the smallest

confidence region (Figure 26), which is marginally larger than the confidence region in Figure 27 corresponding to the reconciliation of six measurements.


Figure 25. Confidence region between kP0 and kT0 when (CM, CT, wP) are reconciled: trended inputs reconciled.

Figure 28. CFMreconciled: contant inputs.

Figure 26. Confidence region between kP0 and kT0 when (Mn, wP, Q, xS) trended inputs are reconciled. Figure 29. CFMestimated: constant inputs.

Figure 27. Confidence region between kP0, and kT0 when (CM, CT, Mn, wP, Q, xS) are reconciled: trended inputs.

Because of the greater information in the measurement and input data, the parameters are estimated with greater confidence than in the case of constant inputs.

The confidence regions are smaller than those obtained for constant inputs. 6.3. EVM Problems. We finally estimated some inputs along with the parameters, giving rise to an EVM problem. Here, the degrees of freedom in NLP 55 are linearly related to the number of finite elements and collocation points. We maintain two collocation points per element and estimate only one input at a time to demonstrate the effect of adding inputs on the reconciled state profiles and parameter estimates. All inputs can be estimated simultaneously. However, NTREST is currently not well suited to this case because of the large increase in the degrees of freedom. We stopped the NLP iterations when the convergence error was less than 10-3, and we added noise from N(0, 0.05 × |unom|), where unom is the value of the inputs reported in Table 1. We then successively estimate CFM, CFT, FF, and QF and hold the remaining inputs constant with 0.5% noise. The results are reported in Figures 28-32, where we plot the reconciled states and inputs. Consider the state profiles in Figure 28 when CFM is estimated. Here, the reconciliation of wP is considerably better than before. Some states are also reconciled closer

Ind. Eng. Chem. Res., Vol. 43, No. 14, 2004 3629 Table 4. Parameter Estimates for Various EVM Problems: Constant Inputs input estimated

parameter estimates {kP0, kT0, kD0}

CFM FF QF CFT

5.05 × 107, 4.51 × 106, 7.87 × 103 5.02 × 107, 4.44 × 106, 7.90 × 103 4.97 × 107, 4.40 × 106, 7.90 × 103 4.95 × 107, 4.38 × 106, 7.87 × 103

Figure 30. FF estimated: constant inputs.

Figure 33. CFMestimated: trended inputs.

Figure 31. QF estimated: constant inputs.

with other inputs estimated. The input profiles show some oscillation toward the end. In addition, the inputs at some of the collocation points were estimated to be at the imposed lower or upper bounds. We also performed an estimation with 60 finite elements, but that did not improve the profiles. In Table 4, we report parameter estimates obtained for the four EVM problems where inputs are constant. To improve the performance of the algorithm and to reduce the disparity in scale of the inputs and parameters, we estimated the logarithm of the Arrhenius factors of the rate constants. Again, the parameter estimates are close to their true values, which is to be expected in this case. We finally performed input estimations for trended inputs, as shown in Figures 33-35. We observe that the inputs seem to have a much better reconciliation, leading to smoother fitted profiles. Table 5 reports the parameter estimates for the EVM parameter estimation with trended inputs. 7. Conclusions

Figure 32. CFTestimated: constant inputs.

to their measurements than before because of the additional parameters in the system. A similar trend in state profiles was also observed for state estimates

We have demonstrated an application of NTREST to parameter estimation of a polymerization reactor model. In particular, we have developed a software interface of NTREST with the DAE2NLP code of Wa¨chter.17 Here, we regard parameters and control variables as independent variables. DAE2NLP uses ADOL-C to generate derivatives of constraints, and we also obtain derivatives of the objective function. In addition, we use basis calculations, matrix inversions, and null-space matrix calculations available from DAE2NLP. The challenges in the estimation problem here are a high-index DAE system, an ill-conditioned discretized DAE, the different information contents of the measurements, and a largescale EVM problem. We reduce the index of the DAE


Acknowledgment Support from the ELKEM Foundation is gratefully acknowledged for this work. The authors also thank Robert Young, Lou Russo, and Robert Fontaine for their help and insights on this problem. Literature Cited (1) Arora, N.; Biegler, L. T. A Trust Region SQP Algorithm for Equality Constrained Parameter Estimation with Simple Parameter Bounds. Comput. Optim. Appl., in press. (2) Albuquerque, J. S.; Biegler, L. T. Data Reconciliation and Gross-Error Detection for Dynamic Systems. AIChE J. 1996, 42, 2841. (3) Narasimhan, S.; Jordache, C. Data Reconciliation and Gross Error Detection; Gulf Publishing Company: Houston, TX, 2000. (4) Crowe, C. M. Observability and Redundancy of Process Data for Steady-State Reconciliation. Chem. Eng. Sci. 1989, 44, 2909. Figure 34.

CFMestimated:

trended inputs.

(5) Zacca, J. J.; Ray, W. H. Modelling of the Liquid-Phase Polymerization of Olefins in Loop Reactors. Chem. Eng. Sci. 1993, 48, 3743-3765. (6) Russo, L. P.; Young, R. E. Moving-Horizon State Estimation Applied to an Industrial Polymerization Process. In Proceedings of the American Control Conference; IEEE Press: Piscataway, NJ, 1999. (7) Bindlish, R. Modeling and Control of Prototypical Polymerization Processes. Ph.D. Thesis, University of Wisconsin, Madison, WI, 1999. (8) Sirohi, A.; Choi, K. Y. On-line parameter estimation in a continuous polymerization process. Ind. Eng. Chem. Res. 1996, 35, 1332-1343. (9) Kiparissides, C.; Seferlis, P.; Mourikas, G.; Morris, A. J. Online optimizing control of molecular weight properties in batch free-radical polymerization reactors. Ind. Eng. Chem. Res. 2002, 41, 6120-6131. (10) Tatiraju, S.; Ray, W. H. Nonlinear state estimation in a polymerization reactor. Ind. Eng. Chem. Res. 1997, 36, 2679-2690. (11) Ramamurthi, Y.; Sistu, P. B.; Bequette, B. W. Controlrelevant dynamic data reconciliation and parameter estimation. Comput. Chem. Eng. 1993, 17, 41-59.

Figure 35. FF estimated: trended inputs.

(12) Cervantes, A. M.; Biegler, L. T. A stable elemental decomposition for dynamic process optimization. Comput. Optim. Appl. 2000, 120, 41-57.

Table 5. Parameter Estimates for Various EVM Problems: Trended Inputs

(13) Biegler, L. T.; Cervantes, A. M.; Wa¨chter, A. Advances in simultaneous strategies for dynamic process optimization. Chem. Eng. Sci. 2002, 57, 575-593.

input estimated

parameter estimates {kP0, kT0, kD0}

CFM QF

4.97 × 107, 4.41 × 106, 7.97 × 103 5.04 × 107, 4.59 × 106, 7.91 × 103

system to make it index 1 and are then able to solve all of the above challenges using NTREST. We generate data by simulating a polymerization reactor. For parameter estimation, we consider several combinations of measured variables whose information content we analyze by constructing joint confidence regions and eigenvalues of the reduced Hessian. Here, we estimate a minimum informative set using four measurements; this is particularly true when the number-average molecular weight (Mn) and the polymer production rate (wP) are measured. Future work will address the development of an extension of NTREST to include a more efficient EVM algorithm that exploits the indefiniteness of the reduced Hessian and an extension of the parameter estimation strategy described in this paper to model discrimination.

(14) Bilardello, P.; Joulia, X.; Lann, J. L.; Delmas, H.; Koehret, B. A General Strategy for Parameter Estimation in DifferentialAlgebraic Systems. Comput. Chem. Eng. 1993, 17, 517. (15) Stewart, W. E.; Henson, T. L.; Box, G. E. P. Model discrimination and criticism with single-response data. AIChE J. 1996, 42, 3055. (16) Stewart, W. E.; Shon, Y.; Box, G. E. P. Discrimination and goodness of fit of multiresponse mechanistic models. AIChE J. 1998, 44, 1404. (17) Wa¨chter, A. An Interior Point Algorithm for Large-Scale Nonlinear Optimization with Applications in Process Engineering. Ph.D. Thesis, Carnegie Mellon University, Pittsburgh, PA, 2001. (18) Mattsson, S.; So¨derlind, G. Index reduction in differentialalgebraic equations using dummy derivatives. SIAM J. Sci. Comput. 1993, 14, 677-692. (19) Hildebrand, B. F. Introduction to Numerical Analysis; Dover Publications: Minneola, NY, 1987. (20) Bader, G.; Ascher, U. A new basis implementation for a mixed order boundary value ODE solver. SIAM J. Sci. Stat. Comput. 1987, 8, 483. (21) Ascher, U.; Spiteri, R.. Collocation software for boundary value differential-algebraic equations. SIAM J. Sci. Comput. 1994, 15, 938-952.

Ind. Eng. Chem. Res., Vol. 43, No. 14, 2004 3631 (22) Griewank, A.; Juedes, D.; Utke, J.; Vogel, O.; Walther, A. ADOL-C: A Package for the Automatic Differentiation of Algorithms Written in C/C++. ACM TOMS 1996, 22, 131. (23) Bard, Y. Nonlinear Parameter Estimation; Academic Press: New York, 1974. (24) Bates, D. M.; Watts, D. G. Nonlinear Regression Analysis and Its Applications; Wiley Series in Probability and Mathematical Statistics; John Wiley and Sons: New York, 1988.

(25) Biegler, L. T.; Nocedal, J.; Schmid, C. A Reduced Hessian Method for Large-Scale Constrained Optimization. SIAM J. Optim. 1995, 5, 314.

Received for review March 13, 2003 Revised manuscript received September 24, 2003 Accepted September 24, 2003 IE030237E

Parameter Estimation for a Polymerization Reactor ... - ACS Publications

Recommend Documents