Abstract

PROCESS SYSTEMS ENGINEERING

A software for parameter estimation in dynamic models

M. Yuceer^I,* * To whom correspondence should be addressed ; I. Atasoy^II; R. Berber^II

^IInonu University Faculty of Engineering Department of Chemical Engineering, Phone: +(90) 422 3410010, Ex: 4753, Fax: +(90) 422 3410046, 44280, Malatya, Turkey. E-mail: myuceer@inonu.edu.tr

^IIAnkara University Faculty of Engineering, Department of Chemical Engineering, 06100, Ankara, Turkey. E-mail: atasoy@eng.ankara.edu.tr, E-mail: berber@eng.ankara.edu.tr

ABSTRACT

A common problem in dynamic systems is to determine parameters in an equation used to represent experimental data. The goal is to determine the values of model parameters that provide the best fit to measured data, generally based on some type of least squares or maximum likelihood criterion. In the most general case, this requires the solution of a nonlinear and frequently non-convex optimization problem. Some of the available software lack in generality, while others do not provide ease of use. A user-interactive parameter estimation software was needed for identifying kinetic parameters. In this work we developed an integration based optimization approach to provide a solution to such problems. For easy implementation of the technique, a parameter estimation software (PARES) has been developed in MATLAB environment. When tested with extensive example problems from literature, the suggested approach is proven to provide good agreement between predicted and observed data within relatively less computing time and iterations.

Keywords: Parameter estimation; Dynamic simulation.

INTRODUCTION

Parameter estimation is a common problem in many areas of process modeling, both in 'on-line' applications such as real time optimization and in 'off-line' applications such as the modeling of reaction kinetics and phase equilibrium. The goal is to determine the values of model parameters that provide the best fit to measured data, generally based on some type of least squares or maximum likelihood criterion. The estimation of parameters in kinetic expressions from time series data is essential for the design, optimization, and control of many chemical systems. The models that describe the kinetics take the form of a set of differential algebraic equations. The statistics and formulation of this parameter estimation problem are well studied (Bard, 1974). In the most general case, this requires the solution of a nonlinear and frequently non-convex optimization problem, which can be approached by dynamic programming (with substantial programming effort and computer time), Pontryagin's maximum principle (that requires solution of adjoint vector whose initial values creating extra problems) or nonlinear programming techniques.

Nonlinear programming may be pursued through discretization of all variables by finite difference approximations or orthogonal collocation. However, inclusion of the collocation variables immensely increases the size of the optimization domain and the degree of polynomial used has a significant impact on the error. The other option is to use the integration approach without transforming the differential equations into a fully algebraic NLP. Pontryagin's two theorems provide continuity and differentiability of the input-output (parameters - state variables) map with respect to the parameters. This method requires calculation of the sensitivities by integrating an extra set of differential equations at the expense of computation time. Esposito and Floudas (2000) provided an underestimating formulation for this method. In their work, two deterministic global optimization methods are introduced for the parameter estimation of models that involve differential algebraic systems. Parameter estimation in dynamic models has been studied and discussed by several workers (Tjoa and Biegler,1991, Albuquerque and Biegler, 1996; Arora and Biegler, 2001, Neumann et al., 2007)

Vassiliadis et al.'s (1994) control vector parameterization method in the same direction, although carries out the optimization in the space of decision variables, uses Lagrange polynomials for expressing control variables and again converts the problem into a finite dimensional NLP with added complexity in the algorithm.

From the perspective of existing software for dynamic model calibration, there are a number of programs available commercially, some of them being accessible via Internet. However, most of them cannot fit some needs of the studied system satisfactorily even though they excel in other fields. For example the OPTKIN software supports large scale sets of equations with many parameters, but little customizability of these equations is allowed (Huybrechts and Van Assche, 1998). Therefore it is quite an efficient tool for the treatment of large sets of first order reactions (e.g. various radical mechanisms) but it does not allow application on some more complicated models. There are also limited possibilities for more detailed analysis of reliability and significance of parameters. The only generally applicable types of programs are the academic and semi-academic ones (Stewart et. al., 1992) but they are supplied as source code and provide poor user interface and require programming skills. Although ERA software package developed by Zamostny and Belohlav (1999) is a useful regression analysis tool, its input data matrix is limited to 20 independent variables and 20 responses, with up to 256 experimental points in each response and the number of model parameters is restricted to 15. Therefore, the floor seems to be open for further developments.

In this work, we have developed an integration based efficient algorithm, and a software for implementation, for parameter estimation in dynamic models without requiring calculation of sensitivity equations. The results obtained for some problems from the literature are compared, and its effectiveness is demonstrated.

OPTIMIZATION OVERVIEW

Optimization techniques are used to find a set of design parameters, x= {x₁, x₂, x₃, .., x_n}, that can in some way be defined as optimal. The objective function, f(x), to be minimized or maximized, may be subject to constraints in the form of equality constraints, inequality constraints and parameter bounds.

General problem may be defined as

subject to the constraints

where f(x), g_j(x) and h_j(x) are scalar functions of the real column vector x. m and r are the number of inequality and equality constraints. g_j(x) and h_j(x) return the values of the equality and inequality constraints evaluated at x, l and u denoting lower and upper parameter bounds. The optimum vector x that solves problem (Eq.1) is denoted by x^* with corresponding optimum function value f(x^*).

Dynamic optimization problem is generally defined as Eq. 3.

x: state variables vector, u: control (optimization) parameters, f: vector function, ψ, h, and g arbitrary functions, p: model parameters, t: time

The fact that dynamic optimization problems are inherently difficult to solve in the first place makes the problem even more challenging. In the solution of dynamic optimization problems the following methods can be used:

▪ Pontryagin's Maximum Principle

▪ Dynamic programming

▪ Converting into a non-linear programming problem (NLP).

The use of maximum principle, which converts the problem into a two-point boundary value problem, poses difficulties in terms of satisfying the constraints represented by the differential equations. Solution by dynamic programming is somewhat problem specific and requires extensive programming effort and time.

Although a wide spectrum of methods exists for unconstrained optimization, methods can be broadly classified whether the derivative information is used or not. Search methods that use only function evaluations (e.g. Nelder-Mead simplex search) are most suitable for problems that are very nonlinear or have a number of discontinuities. Gradient methods are generally more efficient when the function to be minimized is continuous in its first derivative. Higher order methods, such as Newton's method, are only really suitable when the second order information is readily and easily calculated, because calculation of second order information, using numerical differentiation, is computationally expensive (The Mathworks, 2003).

Constrained Optimization

In constrained optimization, the general aim is to transform the problem into an easier subproblem that can then be solved and used as the basis of an iterative process. A characteristic of a large class of early methods is the translation of the constrained problem to a basic unconstrained problem by using a penalty function for constraints that are near or beyond the constraint boundary. In this way the constrained problem is solved using a sequence of parameterized unconstrained optimizations, which in the limit of the sequence converge to the constrained problem. These methods are now considered relatively inefficient and have been replaced by methods that have focused on the solution of the Kuhn-Tucker (KT) equations.

SQP methods represent the state of the art in nonlinear programming methods. At each major iteration, an approximation is made to the Hessian of the Lagrangian function using a Quasi-Newton updating method. This is then used to generate a QP subproblem whose solution is used to form a search direction for a line search procedure. SQP has been extensively described and discussed by several workers (Powell, 1978; Fukushima, 1986, Fletcher, 1987; Boggs and Tolle, 1995), often for large-scale problems, due to their efficiency and robustness. Such a view was based on strong convergence properties of such algorithms, and reinforced in the comparative testing experiments of Hock and Schittkowski (1981).

PROBLEM FORMULATION AND SOLUTION ALGORITHM

The mathematical formulation of chemical reaction mechanisms, taken as exemplary systems to demonstrate the optimization algorithm, is given by a coupled system of stiff nonlinear differential equations

where y is the state vector of the system, p is the model parameters.

In the present work we propose an approach based on converting the dynamic optimization problem into a finite dimensional nonlinear program through discretization of decision variables, like in the method of 'control vector parameterization'. In this case, it is necessary to discretize only the control variables. For given control variables and values of the other variables, it becomes possible to integrate the underlying differential algebraic equation (DAE) system using standard integration algorithms so as to evaluate the objective function and other constraints that have to be satisfied by the solution.

The proposed optimization technique brings some advantages such as ease of implementation, and is effective in the sense that it eliminates any need to evaluate the sensitivity functions, which is an essential part of the computation in optimization strategies based on NLP method.

A control vector parameterization like approach for integration has been suggested such that the search space is discretized in time. By assuming a piecewise constant control, we have been able to integrate the model equations so that the objective function can be evaluated. Running a numerical constrained/unconstrained optimization technique on top of this space allowed us to obtain the values of decision variables, i.e. parameters in the kinetic models, after a reasonable set of iterations to satisfy a stopping criterion. Since no approximation polynomials or sensitivity functions were required, the calculations were substantially simplified.

Figure 1 shows the logic flow diagram for the parameter estimation algorithm developed, which has been extensively tested before with a number of different problems (Agun, 2002). Employing this approach, 33 parameters of water quality model have been identified for a river (Yuceer et al., 2007).

The algorithm offers the possibility of employing different numerical optimization routines with ease to estimate the updated p in order to satisfy the particular needs of the model employed.

Model parameters were estimated by Quasi-Newton (QN), Nelder-Mead Simplex (NMS), Gauss-Newton (GN), Levenberg-Marquardt (LM), Sequential Quadratic Programming (SQP) algorithms by minimizing the objective function, which is the sum of squares of errors between the predicted and measured values for all of the state variables for a dynamic run as follows;

where x: computed value, x₀ : observed value, n: total number of state variables and m: total number of observations. MATLAB^® 6.5 software package and a 2.4 GHz Pentium 4 processor having 512 MB RAM were used for all computations.

PARES SOFTWARE

Software PARES (PARameter EStimation), coded in MATLAB^TM 6.5 has been developed to implement the suggested parameter estimation technique. PARES is an interactive software system to identify parameters in differential algebraic equation system models. The program has ability to make parameter estimation with different optimization methods. For unconstraint optimization problems, Gauss-Newton, Nelder-Mead Simplex, Levenberg-Marquardt or Quasi-Newton methods can be optionally used whereas sequential quadratic programming needs to be used for constraint optimization problems.

The program requires the input data in six main steps described as follows:

1. Model equations (user can enter any model, on the window provided as an M-file)

2. Number of model parameters

3. Initial, lower and upper bounds of the parameters

4. Starting, final and sampling times for experimental data

5. Number of state variables

6. Experimental measurements of state variables

The software, whose graphical user interface is shown in Figure 2, has also the ability to choose different objective functions (i.e. least square error, absolute error or standardized absolute error). It shows optimum model parameters, CPU time, objective function and the fit between experimental and predicted values of each state variable graphically. Contrary to the ERA software package, there is no limit on the number of parameters and state variables.

EXAMPLE PROBLEMS AND RESULTS

Although many linear and nonlinear models with varying degree of kinetic complexity have been tested, three examples are included here briefly to demonstrate the validity of the obtained results with literature for comparison reasons.

Example 1 represents a first-order reversible chain reaction in series, and was presented by Esposito and Floudas (2000). The reaction kinetic is expressed as,

All of the components are supposed to be measured, and therefore their concentrations are included in the model used. The dynamic model is written as

where the state vector, x, is defined as [A, B, C], and parameter vector, p, is defined as [k₁, k₂, k₃, k₄].

Figure 3 depicted that the comparison of experimental and observed data using SQP for Example 1. The results were summarized in Table 1.

Example 2 tackles catalytic hydrogenation of cinnamaldehyde, which was previously investigated by Zamostny and Belohlav (1999). Reactions occurring in the system are indicated in the scheme as follows. The experiment was carried out in an isothermal, isobaric, stirred semi-batch reactor, the mass balance of which is given by the following equations.

where r : reaction rate, k : rate constant, t : time, and θ_i : catalyst surface coverage for i-th compound. The parameter estimation problem was solved using different optimization routines.

Figure 4 shows that the predicted data fits, in a good agreement, with the experimental measurements. By increasing the complexity of the problem, SQP method becomes more effective than compared to those other methods. The results were summarized in Table 2.

Example 3 The kinetic is a first-order irreversible chain reaction as presented by Esposito and Floudas (2000).

Only the concentrations of components A and B were measured; therefore, component C did not appear in the model used for estimation. The differential equation model takes the form

where the state vector, x, is defined as [A, B], and parameter vector, p, is defined as [k₁, k₂]. The data used in the study was generated with the p values of [5, 1].

The parameter estimation problem was solved using different optimization methods. Figure 5 shows that the observed data fits the experimental measurements very well. The results were summarized in Table 3. According to these results, the number of iteration and CPU time of SQP method are much better than the other methods.

CONCLUSIONS

As reflected by the results, solution to such difficult dynamic optimization problems for parameter estimation seems to be effectively achieved by suggested approach. Particularly, Table 1 reveals that PARES software, which is capable of accomplishing the task with much less iterations and CPU time than previously published techniques, and furthermore the SQP method appears to be favourable to other numerical optimization routines.

(Received: April 1, 2007; Accepted: May 5, 2008)

Publication Dates

History