- Email: [email protected]

1091

Developing Phenomena Models from Experimental Data Niels Rode Kristensen", Henrik Madsen^ and Sten Bay J0rgensen^ "Department of Chemical Engineering, ^Informatics and Mathematical Modelling, Technical University of Denmark, DK-2800 Lyngby, Denmark

Abstract A systematic approach for developing phenomena models from experimental data is presented. The approach is based on integrated application of stochastic differential equation (SDE) modelling and multivariate nonparametric regression, and it is shown how these techniques can be used to uncover unknown functionality behind various phenomena in first engineering principles models using experimental data. The proposed modelling approach has significant application potential, e.g. for determining unknown reaction kinetics in both chemical and biological processes. To illustrate the performance of the approach, a case study is presented, which shows how an appropriate phenomena model for the growth rate of biomass in a fed-batch bioreactor can be inferred from data.

1. Introduction For most chemical and biological processes first principles engineering methods can be applied to formulate balance equations that essentially provide the skeleton of an ordinary differential equation (ODE) model for such a process. What often remains to be determined, however, are functional relations in the constitutive equations for phenomena such as reaction rates and heat and mass transfer rates. These phenomena models are often difficult to determine due to the fact that finding a parametric expression with an appropriate structure to match the available experimental data is essentially a trial-and-error procedure with limited guidance and therefore potentially time-consuming. In the present paper a more systematic procedure is proposed. The key idea of this procedure is to exploit the close connection between ODE models and SDE models to develop a methodology for determining the proper structure of the functional relations directly from the experimental data. The new procedure more specifically allows important trends and dependencies to be visually determined without making any prior assumptions and in turn allows appropriate parametric expressions to be inferred. The proposed procedure is a tailored application of the grey-box modelling approach to process model development proposed by Kristensen et al. (2002b), within which specific model deficiencies can be pinpointed and their structural origin uncovered to improve the model. The remainder of the paper is organized as follows: In Section 2 the details of the proposed procedure are outlined; in Section 3 a case study illustrating the performance of the procedure is presented and in Section 4 the conclusions of the paper are given.

1092

ODE Oiodel formutalKW

^ ^

SDEitJOdeJ

^

State esumalton

^ ^

Ncmparame^dc mcRfeB«i&

A

t

f

First engineering principles

Experimental data

Estimate of functional relation

Figure 1. The proposed procedure for developing phenomena models. The boxes in grey illustrate tasks and the boxes in white illustrate inputs and outputs.

2. Methodology The proposed procedure is shown in Figure 1 and consists offivebasic steps. First a standard ODE model is derived fromfirstengineering principles and the constitutive equations containing unknown functional relations are identified. The ODE model is then translated into a stochastic state space model consisting of a set of SDE's describing the dynamics of the system in continuous time and a set of discrete time measurement equations signifying how the available experimental data was obtained. A major difference between ODE's and SDE's is the inclusion of a stochastic term in the latter, which allows uncertainty to be accomodated, and which, if the constitutive equations of interest are reformulated as additional state equations, allows estimates of the corresponding state variables to be computed from the experimental data. The specific approach used for this purpose involves parameter estimation and subsequent state estimation by means of methods based on the extended Kalman filter (EKF). By subsequently applying methods for multivariate nonparametric regression to appropriate subsets of the state estimates, visual determination of important trends and dependencies is facilitated, in tum allowing appropriate parametric expressions for the unknown functional relations in the constitutive equations to be inferred. More details on the individual steps of the proposed procedure are given in the following. 2.1. ODE model formulation In the first step of the procedure, a standard ODE model is derived and the constitutive equations containing unknown functional relations are identified. Deriving an ODE model from first engineering principles is a standard discipline for most chemical and process systems engineers and in the general case gives rise to a model of the following type: dxt

"df

f(xt,ut,rt,t,0)

(1)

where t G Eistime,a;t G W^ is a vector ofbalanced quantities or state variables, u^ G W^ is a vector of input variables and ^ G M^ is a vector of possibly unknown parameters, and where /(•) G E"* is a nonlinear function. In addition to (1) a number of constitutive equations for various phenomena are often needed, i.e. equations of the following type:

n = ^{xt,ut,6)

(2)

1093 where n is a given phenomenon and (/?(•)€ M is the nonlinear function of the state and input variables needed to describe it. This function is, however, often unknown and must therefore somehow be determined from experimental data. In the context of the systematic procedure proposed in the present paper, the first step towards determining the proper structure of (^(•) is to assume that this function and hence rt is constant. 2.2. SDE model formulation In the second step of the procedure the ODE model is translated into a stochastic state space model with rt as an additional state variable. This is straightforward, as it can simply be done by replacing the ODE's with SDE's and adding a set of discrete time measurement equations, which yields a model of the following type: dxl = f*{xt,Ut,t,e)dt-\-a*{ut,t,e)duj* Vk =h{xl,Uk,tk,9)-\-ek

(3) (4)

where t G E is time, x^ = [xj r^]^ G E"^"*"^ is a vector of state variables, Ut € M"* is a vector of input variables, y^ € EMS a vector of output variables, ^ G E^ is a vector of possibly unknown parameters, /*(•) G E ^ + \ or*(-) G E^^+^^^^^^^+i) and /i(-) G Mf are nonlinear functions, {a;^} is an (n -f 1)-dimensional standard Wiener process and {e^} is an /-dimensional white noise process with e^ G AT (0, S(uk^tk^ 9)). The first term on the right-hand side of the SDE's in (3) is called the drift term and is a deterministic term, which can be derived from the term on the right-hand side of (1) as follows: /V^* ., + a\ J {xt,ut,t,e) =

1

(f{xt,uurt,t,e)\ 1

(5)

where the zero is due to the assumption of constant r^. The second term on the right-hand side of the SDE's in (3) is called the diffusion term. This is a stochastic term included to accommodate uncertainty due to e.g. approximation errors or unmodelled phenomena and is therefore the key to subsequently determining the proper structure of (^(•). A more detailed account of the theory and application of SDE's is given by 0ksendal (1998). 2.3. Parameter estimation In the third step of the proposed procedure the unknown parameters of the model in (3)(4) are estimated from available experimental data, i.e. data in the form of a sequence of measurements yQ,y^,...,yj^,...,y^. The solution to (3) is a Markov process, and an estimation scheme based on probabilistic methods, e.g. maximum likelihood (ML) or maximum a posteriori (MAP), can therefore be applied. A detailed account of one such scheme, which is based on the EKF, is given by Kristensen et al. (2002a). 2.4. State estimation In the fourth step of the procedure, state estimates are computed to facilitate determination of the proper structure of (^(•) by means of subsequent multivariate nonparametric regression. Using the model in (3)-(4) and the parameter estimates obtained in the previous step, state estimates xl^j^, A: = 0 , . . . , iV, can be obtained by applying the EKF once again using the same experimental data. In particular, since vt is included as an additional

1094 state variable in this model, estimates r^jj^, A: = 0 , . . . , AT, can be obtained, which in turn facilitates application of multivariate nonparametric regression to provide estimates of possible functional relations between rt and the original state and input variables. 2.5. Nonparametric modelling In the fifth step of the procedure the state estimates computed in the previous step are used to determine the proper structure of (p{-) by means of multivariate nonparametric regression. Several such techniques are available, but in the context of the proposed procedure, additive models (Hastie and Tibshirani, 1990) are preferred, because fitting such models circumvents the curse of dimensionality, which tends to render nonparametric regression infeasible in higher dimensions, and because results obtained with such models are particularly easy to visualize, which is important. Additive models are nonparametric extensions of linear regression models and are fitted by using a training data set of observations of several predictor variables X i , . . . , Xn and a single response variable Y to compute a smoothed estimate of the response variable for a given set of values of the predictor variables. This is done by assuming that the contributions from each of the predictor variables are additive and can be fitted nonparametrically using the backfitting algorithm (Hastie and Tibshirani, 1990). The assumption of additive contributions does not necessarily limit the ability of additive models to reveal non-additive functional relations involving more than one predictor variable, since, by proper processing of the training data set, functions of more than one predictor variable, e.g. X1X2, can be included as predictor variables as well (Hastie and Tibshirani, 1990). Using additive models, the variation in ffc|jfc, A: = 0 , . . . , A/^ can be decomposed into the variation that can be attributed to each of the original state and input variables and the result can be visualized by means of partial dependence plots with associated bootstrap confidence intervals (Hastie et ai, 2(X)1). In this manner, it may be possible to reveal the true structure of (p{-) and subsequently determine an appropriate parametric expression for the revealed functional relation. Remark. Once an appropriate parametric expression for the unknown functional relation has been determined, the parameters of this expression should be estimated from the experimental data and the quality of the resulting model should subsequently be evaluated by means of cross-validation. A discussion of methods for evaluating the quality of a model with respect to its intended application is given by Kristensen et al. (2002b). Remark. A key advantage of the proposed procedure is that functional relations involving unmeasured variables can easily be determined as well if certain observability conditions are fulfilled, e.g. functional relations between reaction rates, which can seldom be measured directly, and concentrations of various species, which may also be unmeasurable.

3. Case Study: Determining the Growth Rate in a Fed-batch Bioreactor To illustrate the performance of the proposed procedure, a simple simulation example is considered in the following. The process considered is a fed-batch bioreactor, where the true model used to simulate the process is given as follows:

(6)

1095 where X is the biomass concentration, S is the substrate concentration, V is the volume of the reactor, F is the feed flow rate, Y is the yield coefficient of biomass, 5 ^ is the feed concentration of substrate, and /x(5) is the biomass growth rate, which is characterized by Monod kinetics and substrate inhibition, i.e.: l^[S) = /in

(7)

'K2S'^-^S-\-Ki

where /Xmax, Ki and K2 are kinetic parameters. Simulated data sets from two batch runs are generated by perturbing the feed flow rate along a pre-determined trajectory and subsequently adding Gaussian measurement noise to the appropriate variables (see below). Using these data sets and starting from preliminar balance equations, where the biomass growth rate is assumed to be unknown, it is illustrated how the proposed procedure can be used to visually determine the proper structure of /i(5). In the context of the first step of the proposed procedure, an ODE model corresponding to (6) has thus been formulated and the constitutive equation for /i(5) has been identified as containing an unknown functional relation. The first step towards determining this relation then is to assume that /i(5) is constant, i.e. ^{S) = /JL, and translate the ODE model into a stochastic state space model with /i as an additional state variable, i.e.:

dt +

+ eik,

•(Til

0

0

a22

0 .0

0 0

ekeN{0,S)).

s=

0 0

0

0" 0 da,; 0

(8)

•5ii

0

0 . 0

522

0 •

0

0

5'33.

(9)

As the next step, the unknown parameters of this model are estimated using the EKFbased estimation scheme presented by Kristensen et al. (2002a) and the data sets mentioned above. Using the resulting model and the same data sets, state estimates X^^, Sk\k^ yk\k^ P'kiky A; = 0 , . . . , A^, are then obtained by applying the EKF once again and an additive model is fitted to reveal the true structure of the function describing /x by means of estimates of functional relations between fi and the state and input variables. It is reasonable to assume that /i does not depend on V and F, so only functional relations between fik^k ^^^ ^k\k and 5^]^; are estimated, giving the results shown in Figure 2. These plots indicate that /lk\k does not depend on XJ^JJ^, but is highly dependent on Sk\ky which in turn suggests to replace the assumption of constant /x with an assumption of /x being a function of S that complies with the revealed functional relation. To a skilled biotechnologist this relation is a clear indication of Monod kinetics and substrate inhibition and immediately suggests to replace /x with the true function fj,{S) in (7). Remark. The case study presented here illustrates how the proposed procedure can be used to determine the proper structure of an unknown functional relation direcdy from experimental data by providing information about which variables contribute to the observed variation and how these contributions may be characterized parametrically. Once a parametric expression has been determined, the parameters can be estimated and the

1096 0.25

0.2 0.15 ai

0.05

..... :^r.^ .},'.•••••''•

B-

•. .jji>-H^^ ^ • r

0 •y /

•

•'•'

"i

.•••^

•. •

-0.05 -0.1 -0.1S

-02

1.5

2

2.6

3

(a)Afc|fc vs.Xfciib

3.5

(b)Afc|fc vs. 5;fe|fc

Figure 2. Partial dependence plots of [jLk\k vs. X^^k and Sk\k (solid lines: Estimates; dotted lines: 95% bootstrap confidence intervals).

quality of the resulting model evaluated by means of cross-validation. The proposed procedure thus facilitates systematic development of phenomena models from data.

4. Conclusion A systematic approach for developing phenomena models for chemical and biological processes from experimental data has been presented. The approach is based on integrated application of stochastic differential equation (SDE) modelling and multivariate nonparametric regression and can be used to uncover unknown functionality behind various phenomena in first engineering principles models using data. The proposed modelling approach has significant application potential, e.g. for determining unknown reaction kinetics in both chemical and biological processes, and a key advantage in this regard is that functional relations involving unmeasured variables can easily be determined as well if certain observability conditions are fulfilled.

5. References Hastie, TJ. and Tibshirani, RJ. (1990). Generalized Additive Models. Chapman & Hall, London, England. Hastie, TJ.; Tibshirani, RJ. and Friedman, J. (2001). The Elements of Statistical Learning - Data Mining, Inference and Prediction. Springer-Verlag, New York, USA. Kristensen, NR.; Madsen, H. and J0rgensen, SB. (2002a). Parameter Estimation in Stochastic Grey-Box Models. Submitted for publication. Kristensen, NR.; Madsen, H. and J0rgensen, SB. (2002b). A Method for Systematic Improvement of Stochastic Grey-Box Models. Submitted for publication. 0ksendal, B. (1998). Stochastic Differential Equations - An Introduction with Applications. Springer-Verlag, Berlin, Germany, fifth edition.

Copyright © 2021 COEK.INFO. All rights reserved.