- Email: [email protected]
Statistical tools for optimal dynamic model building
ELSEVIER
Computers and Chemical Engineering 24 (2000) 1261-1267
Computers &Chemical Engineering www.elsevier.com/locate/compchemeng
Statistical tools for optimal dynamic model building S.P. Asprey *, S. Macchietto Centre for Process Systems Engineering, Imperial College of Science, Technology and Medicine, Prince Consort Road, Londoln SW 7 2BY, UK
Abstract
A general, systematic procedure is presented to support the development and statistical verification of dynamic process models. Within this procedure, methods are presented to address several key aspects, such as structural identifiability and distinguishability testing, as well as optimal design of dynamic experiments for both model discrimination and improving parameter precision. A novel optimisation-based approach is introduced for testing of model structural identifiability and distinguishability, involving semi-infinite programming and max-rain problems. The design of dynamic experiments is cast as an optimal control problem within a framework that enables the calculation of optimal sampling points, experiment duration, fixed and variable external control profiles, and initial conditions of a dynamic experiment subject to general constraints on inputs and outputs. Within this framework, methods are presented to provide experiment design robustness, accounting for parameter uncertainty. The procedure is demonstrated through a practical biotechnology example. © 2000 Elsevier Science Ltd. All rights reserved. Keywords: Optimal experiment design; Dynamic modelling; Parameter estimation
1. Introduction
Detailed mathematical models are increasingly used by engineers in academia and industry to gain competitive advantage through such applications as modelbased process design, operations, and control. Thus, building and verifying high quality steady-state or dynamic, single or multiresponse mechanistic models of processing systems are key activities in process engineering. In this work, we present a general, systematic procedure to support the development and statistical verification of dynamic process models. The procedure comprises components that address several key aspects, such as structural identifiability and distinguishability testing, and optimal design of dynamic experiments for model discrimination and improving parameter precision. Experiment design is concerned with the following question: How does one adjust the time-varying controis, initial conditions or length of experiment for the experiment to generate the maximum amount of infor-
* Corresponding author. Tel.: + 44-20-75946604; fax: + 44-2075946606. E-mail address: [email protected] (S.P. Asprey)
mation for the purpose of determining the correct model from a set of candidates, and in turn, for estimating the parameters with greatest precision? For this, the problem is cast as an optimal control problem. Within this framework, one can calculate optimal sampiing points, experiment duration, fixed and variable external controls with input constraints, and initial conditions. To mathematically represent time-varying external controls to the process, we use the control vector parameterisation (CVP) technique, with either piecewise constant, or higher-order functions with 0th or higher-order continuity (Vassiliadis, Sargent & Pantelides, 1994). Here, we first present the overall procedure, followed by a brief discussion of current and new methods for each key component. Potential applications of these techniques in all engineering disciplines are abundant, including applications in chemical kinetics and reaction mechanism elucidation, polymer reaction engineering, physical properties estimation, and biochemical engineering. Throughout the paper, an engineering example using a typical semicontinuous bioreactor application is used to demonstrate the overall procedure and the power of the proposed techniques in their ability to reduce the quantity of experimental work required, while increasing the quality of the results.
0098-1354/00/$ - see front matter © 2000 Elsevier Science Ltd.'All rights reserved. PII: S0098-1354(00)00328-8
S.P. Asprey, S. Macchietto / Computers and Chemical Engineering 24 (2000) 1261-1267
1262
2. The modelling process When building models whose use is to explain an observed phenomenon, one uses a priori knowledge such as physical, chemical or biological 'laws' to propose (conceivably several) possible models. In each case, these laws dictate the model structure, and we may wish to know whether one or more such structures are adequate for the problem at hand. These models contain parameters that may have physical meaning, and we may wish to know if it is at all possible to determine their values, and, if so, to do so with maximum precision. In what follows, we therefore consider general deterministic models in the form of a set of (possibly mixed) differential and algebraic equations: f(~(t), x(t), u(t), w, 0, t) = 0 y(t) = g(x(t))
(1)
where x(t) is an ns-dimensional vector of time dependent state variables, u(t) is an no-dimensional vector of time-varying controls or inputs to the process, w is an n,.-dimensional vector of constant controls, 0 is a P-dimensional vector of model parameters to be determined within a continuous, realisable set 19, and y(t) is an M-dimensional vector of measured response variables that are a function of the state variables, x(t). We first present an overall model building strategy, depicted in Fig. 1. The first step to this strategy involves specifying one or more models to describe the process at hand. In Stage I, the strategy then uses these models to perform preliminary structural identifiability and distinguishability tests before any data are collected. These tests are used to determine whether or not the parameters in the mathematical formulation of the models can be uniquely identified, and whether or not the model structures can be distinguished from one another. In Stage II, the strategy involves designing experiments for model discrimination, followed by ..... .r..,.p..~.~.ry...A..,r.].y:!.~ ........ ~..~.fl .~ .,...a~....,r.n..,r. ....... .r.,.~..~.~..r...,~...L,., " ....
if* ©
E
parameter estimation and model adequacy checking to determine whether or not any of the proposed models can be rejected, with the final purpose of selecting the 'best' model from the given set. Once we are left with the 'best' model, we then go on to Stage III to design experiments for improving the precision of the parameters within that model, eventually arriving at a final, statistically-verified model formulation.
3. Preliminary analysis 3.1. Structural identifiability and distinguishability During the preliminary model analysis, we simply wish to know whether, given a candidate set of models and a proposed set of variables to be measured, we can in principle design an experiment such that the parameters of these models can be uniquely identified (identifiability), and whether the different model structures can be distinguished from one another (distinguishability). Past approaches for testing identifiability in nonlinear-in-input dynamic models find their roots in differential geometry, where, in particular, the model is approximated locally by its input-output mapping using functional expansions (Walter, 1987). Coefficients of these functional expansions are then evaluated at parameter values 0 and 0", and used to form systems of nonlinear algebraic equations, whose unique solution 0 = 0 " ensures identifiability. This approach has its advantages and disadvantages; its advantages being (i) it is a definitive test (answers the test with yes/no), and (ii) the methods are not affected by scaling problems. Disadvantages of the methods are (i) they are limited to models of special forms (e.g. nonlinear control-affine models), and (ii) they are limited to models with dimensions (ns+ P)_< 10, very small for modelling realistic process systems. To avoid these disadvantages (in particular disadvantage (ii)), we propose a novel optimisation-based approach, defined in the following (Pantelides, 2000).
Definition: Identifiability. A model giving output trajectory y(0, u(t), t) is globally identifiable if, for any two parameter sets, 0 c O and 0"~O*, a time horizon of interest, t~[0, tA, for all system inputs u(t)eU, and the same initial conditions y(t = 0), the global maximum: Test II:
~z=
max
0~0,0"~0"
(0-0*)zW0(0-0*)
(2a)
subject to: (y(u(t), O) -- y(u(t), O*))rWy(y(u(t), O) -- y(u(t), 0")) Fig. 1. The overall model-building scheme.
x dt < eyV u(t)eU
S.P. Asprey, S. Macchietto / Computers and Chemical Engineering 24 (2000) 1261-1267 12' 1' 0.8,
0.2, I ti=100
I
0 9a
40
60
80
I00
I~
Fig. 2. Structural identifiabilityand distinguishabilityillustration.
f(~(t), y(t), x(t), u(t), w, 0, t) = 0 (the model equations)
0~ < 0 , < 0~v 07L < 0 * < 0 T v
i=1 ..... P (2b)
Definition: Distinguishability. Two models giving output trajectories y(O, u(t), t) and y*(0*, u(t), t), respectively, are distinguishable if for any 0~O, 0*eO*, re[O, t/] for all u(t)eU, and the same initial conditions y(t = 0), the global maximum: Test DI:
fT min
u(t)EU0~0,0"~0" 0
(yo)(t) - y(2)(t)) r
x W(yo)(t ) - y{2)(t)) dt
(3a)
fo)(£(t), y(t), x(t), u(t), w, 0, t) = 0 f(2)(~,(t), y(t), x(t), u(t), w, 0", t) = 0
i=l,...,P
0 * L < 0 * < 0 *t~ i = l . . . . . P
4.1. Model d&crimination The problem of model discrimination occurs when ArM rival models have been proposed to describe a system and it is not certain which of the ArM models is 'best'. For any experiment, it is expected that the 'best' model will provide the most accurate prediction of the observed measurements. Experimental designs optimised for model discrimination are termed T-optimal, and are based on maximising an overall measure of the divergence between the model predictions. Very little past work has been presented for discriminating between rival dynamic models. Espie and Macchietto (1989) were the first to attempt an extension of the Hunter and Reiner (1965) criterion to the dynamic case, where the summation over the data points in the original formulation was replaced by an integral over a pre-specified time interval [to, t/]: ft/~ Um ^^ q~r = arg max ~ Tt,f(g~,O, 0', t) dt (4) r,,,,(~,
subject to:
Of
responses of the two models, over the time horizon of interest, If such a value is sufficiently large, the two models may be deemed distinguishable. Again, returning to our illustration of Fig. 2, we wish to ensure that when the same input u(t) is applied to two different models with parameter vectors 0 and 0", respectively, and the same initial conditions, the two predicted model output trajectories do not overlie. This optimisation approach is not limited to any particular model form or size; however, we are now faced with solving semi-infinite or m a x - m i n dynamic optimisation problems. We must also choose a priori tolerances for each test (E% and ~ , best chosen to be commensurate with the integral o~ the expected variance of experimental measurements). A more complete derivation of the optimisation-based approach to testing identifiability and distinguishability, as well as solution methods will be given in a future paper.
4. Design of dynamic experiments
i=l ..... P
gives @~< %. Here e%and ey are arbitrarily small positive values and W0 and Wy are weighting matrices. In words, @z is the largest distance between two parameter vectors that still give (essentially) the same response trajectories. If such a distance is arbitrarily small, the model can be deemed globally identifiable. In our illustration of Fig. 2, we wish to ensure that when the same input u(t) is applied to the same model, at the same initial conditions, but with different parameter vectors 0 and 0* (such that 0 ~ 0"), the two predicted model output trajectories do not overlie.
@D = max
1263
(3b)
gives @D > ~ o for some ~ o > 0. In words, @D is the largest value of the minimum distance between the
~, ~', t)
=
(~,(e, ~, t)
- ~,,(to, g', t)) ~
× W(~t(q), 0, t) -- ~r(q), 0', t))
(5)
and embedded in an optimal control framework. Here, is the vector of the best available estimates of the model parameters, and (p is the vector of experiment design variables. However, by using a continuous function (Eq. (4)), the meaning of the sampling time of the responses is lost. In our newly proposed formulation, we keep the concept of sampling time by using the Hunter-Reiner formulation as a starting point for the design of dynamic experiments, where now we use the notion of discrete sampling points, tk, k = 1 , . . . , nsp of each of the response measurements:
S.P. Asprey, S. Macchietto / Computers and Chemical Engineering 24 (2000) 1261-1267
1264
arg max
~1[~T =
NM
NM
~,
~
T4r(~O,O, 0', tk)
(6)
~p~O k = l l = l 1'=l--1
We embed this criterion within a dynamic optimisation framework, described in detail in Asprey, Macchietto and Pantelides (2000). The design vector, q~, now comprises the discrete sampling points as well as any time-varying and time-invariant inputs, and variable initial conditions. To account for uncertainty in the point estimates of the parameters, and thus variance in the predictions of the models, new extensions of (6) could be formulated based on the Bayesian posterior density approach used for nonlinear steady-state process models (Stewart, Shon & Box, 1998), or by formulating either an expected value criterion: qJTEV
= arg max q~
f ,~= NM
E
"~ ~
O~O,O'EO' k
NM
1 l=1 l'=1+1
--}
min ~[V(0, ¢p, ws)]. Cp,W
Tu,(~p, 0, 0', tk)
(7) or a worst-case criterion:
tPrwc = arg max ~
min
NM
NM
Y',
~
However, this formulation is restricted to a diagonal variance-covariance matrix of the parameters, an assumption that is almost always false. Again, the continuous function (Eq. (11)) approach does not explicitly allow one to optimise discrete sampling times of the system responses. More recent work has been carried out by KSrkel, Bauer, Bock and Schl6der (1999) on sequential experiment design for parameter precision improvement in models described by DAEs. They take a similar optimal control approach, using piecewise constant or piecewise linear functions for parameterising the time-varying controls, but instead use a metric (@[] = determinant, D-optimal; maximum eigenvalue, E-optimal, etc.) of the parameter variance-covariance matrix as an objective function:
Tt,r(~p, 0, @, tk)
OeO,O'~O' k = 1 1=1 l ' = l + 1
(8)
For selecting optimal sampling times of the system responses, they use binary decision variables, ws, as weights to select from a set of equally spaced points along the time horizon of interest. Much of our work to date is based on a unique criterion put forth by Zullo (1991), who defined a discrete version of the dynamic information matrix for experiment design for the improvement of parameter precision, defined by:
4.2. Improving parameter precision
M
=
The aim of experiment design for improved parameter precision is to decrease the size of the inference regions of each of the parameters in any given model. This is equivalent to making the elements of the parameter variance-covariance matrix 'small'. Thus, we use the marginal posterior density covariance, given by (Box & Lucas, 1959): V(O, q~) :'~
[~ f r
Vl'~'rs~T~', ~'~
]--1
=
[ Srl
Sr2
Qr = ISrl
(10)
~r{~ is the rsth element of the inverse of the estimate of the variance-covariance matrix ~y = cov(y,, y~). Espie and Macchietto (1989) were the first to attempt an extension of these concepts to the dynamic case for systems described by DAEs, formulating the design of experiments as an optimal control problem and suggesting a criterion based on maximising the weighted sum of the squared dynamic sensitivity coefficients: ~ wi/(S,y)2 dt
qJ~ ¢l~ .it 0 j = 1 i ~ 1
""
Srp
] ~
S,2 "'" S~e[ ~
1 S,2 •
80
q~ = arg max
where tp is the vector of experiment design variables, as above. The matrix Qr is the matrix of dynamic sensitivity coefficients of the rth response variable in the model computed at each of the n~e sampling points:
ls=1
8y,(cp,0, t) o=~
(12)
r=ls=l
sampling point sampling point
t~p1 t~p2
sampling point
tsp,sp
(9)
where Sr is the vector of partial derivatives of the rth response variable in the model with respect to the parameters 0 calculated at all experimental points: S,
M
E
(11)
"'" ,
Sr e J
~
(13) Extensions of this criterion to account for the uncertainty in the point estimates of the parameters are made in Asprey et al. (2000), in which the D-optimal criterion is embedded within a robust optimisation framework, using either an expected value approach, based upon a multivariate normal probability distribution function (pdf) describing the prior uncertainty in the parameters, 0 ,,~ N(O, 2~):
q~ev = arg max E {[Mz(0, q~)]} q~
(14)
0~
where: oE { I M / ( o,
,p)I}=
IM/(0, ¢p)
S.P. Asprey, S. Macchietto / Computers and Chemical Engineering 24 (2000) 1261-1267
or a worst-case approach over an admissible domain O: • R = arg max min {[MI(0, tp)[} tpe~,
(15)
5. A process engineering example
Here, an example is presented to illustrate the application of the dynamic model building tools presented in Sections 2-4. Consider a fed-batch reactor (see Fig. 3) in which the fermentation of baker's yeast is carried out. Within this process, there are two time-varying controls (ul, Uz) and two measured concentrations
1265
5.1. Identifiability o f Mz(O, u, t)
Global structural identifiability of MI(O, u, t) can be tested via Test I1 and using solution techniques for such semi-infinite programming problems (see Asprey et al., 2000). The solution, using 0~e[0.01, 0.98] Vi, a five-element piece-wise constant CVP for the process inputs and a time horizon of interest of 40 h, was found to be ~ y = 1.166 x 10 -6 ( < g.a,t --- 1.0 x 10-5; gy :- 4.0 x 10-3), n o t significantly different from zero and thus indicating the global identifiability of 34/(0, u, t). 5.2. Distinguishability o f MI(O, u, t) and MIr(O, u, t)
( X l , X2).
We propose two models to describe the system. Two differential equations, common to both models, describe the component balances within the reactor: dXl dt = (r - Ua - 04)xl dx--3 = - rx--2+ ul(u2 - x2)
dt
(16)
03
The first model assumes Monod-type kinetics for substrate consumption, giving a third model equation: MI(0, u, t):
01X2
r=
02 + x2
(17)
The second model assumes Canoid-type kinetics and thus, instead of Eq. (17), gives the model form: MII(O, u, t):
01x2
r=
(18)
0 2 X 1 -4- X 2
The experimental conditions that characterise a particular experiment are: (i) the initial biomass concentration (or inoculation), x °, with range 1-10 g 1-~; (ii) the dilution factor, ua, with range 0.05-0.20 h-1; (iii) the substrate concentration in the feed, u2, with range 5-35 g 1-1. The initial substrate concentration, x2(0), is always 0.01 g 1-~ and cannot be manipulated for experiment design purposes. Both x~ and x2 can be measured during the experiment (y = [xb x2]r).
nip~uts
a~-v~n¢"'~,_
Global structural distinguishability of Mz(0, u, t) and Mn(O, u, t) can be tested by reformulating the m a x -
min Test D1 problem to a semi-infinite programming problem and using the techniques mentioned above. The resulting solution of the distinguishability test using the same conditions and CVP as above was found to be ~ D = 3.392 x 10 -1, giving a maximum larger than e~o.= 4.0 × 10 - 3 and thus confirming the distinguishabdity between MI(0, u, t) and MzI(O, u, t). 5.3. Designing a dynamic experiment for model discrimination
Using models 3/1i(0, u, t) and MI(0, u, t), we begin by assuming availability of initial parameter estimates 01 = [0.30, 0.25, 0.56, 0.02] r and On = [0.30, 0.03, 0.55, 0.03] r, respectively. A dynamic experiment design was then carried out for model discrimination using Eq. (6) and piece-wise constant control inputs, with a maximum experiment duration of 40 h. An a priori selection of the number of sampling points, along with the number of switching times/constant levels of the input profiles is required. We arbitrarily chose ten sampling times (nsp = 10) and five time intervals for the control vector parameterisation (Gwi = 5, i = l, " ' " ' n, where, for the • bloreactor models nu 2). We start the optimisation with an experiment lasting I0 h, with hourly sampling points and constant control inputs at 0.2 and 15.0), giving an initial divergence of 1.98. We impose the constraints (cf. Asprey et al., 2000): 1.0
_
5.0
Fig. 3. A fed-batch reactor.
i =
1. . . . ,nsv
i = 1 . . . . . G;
j = 1. . . . . Gw.~ The optimal values for the design vector maximising criterion (6) are presented diagrammatically in Figs. 4 and 5, giving an optimal divergence of 1.19 x 105. The model predictions from the optimal design are shown in Fig. 4. The optimal experiment lasts 40 h and vertical arrows along the abscissa show the optimal sampling points of the two response variables (both measured at the same sampling points).
S.P. Asprey, S. Macchietto / Computers and Chemical Engineering 24 (2000) 1261-1267
1266
By simulating additional 'data' at these experimental conditions, with a homoscedastic normal error variance of 0.04 for both responses, subsequent parameter estimation gives the results in Table 1.By inspection of the multiresponse Z 2 lack-of-fit analyses shown in Table 1, we can conclude that model MI(0, u, t) is better suited to describe the observed experimental data, although, statistically, both models describe the observed data very well.
111.0
~
10•0
'
14,0
'
12•0
'
0.0
/'0.0
YOI(21- . , . . ~ . 5.0
,.2.0
.
. 10.0
.
I~
tt
. 15.0
20.0
25.0
Time
30.0
35.0
45.0
40.0
[11]
5.4. Designing a dynamic experiment for improving parameter precision
Fig. 4. The model predictions for model discrimination.' 0.3.
36.0,
£0.2-
g 3S.0,
LU
li.o 0.1" ! 0 . 1 - ~ 0.0 10.0
0.0
~ 31.0'
,
_, : _,
20.0
30.0
40.0
50.0
300 0.0
10.0
TmN
- - 20.0 30.0 llmet~l
40.0
50,0
Fig. 5. The optimal input profiles for model discrimination.
Table 1 Parameter estimates for MI(0, u, t) and Mn(0, u, t) after one designed experiment Parameter
Estimate
Confidence interval
0= 02 03
Mz 0.325 0.267 0.559 0.055
Mu 0.303 0.011 0.549 0.049
Mz 40.012 +0.075 __40.018 + 0.006
0.58 1.96
Z2ef= 26.3 Z2ef= 26.3
04
LOF (M/) LOF (Mzz)
MI! _+0.011 +0.009 _+0.023 _+0.011
We next turn to improving the precision in the parameter estimates in MI(O, u, t). Using the estimates shown in Table 1, an experiment design for improving parameter precision using the E-optimality criterion with Eq. (12) was performed. Using the same initial values, constraints, and representation of the time-varying inputs as in the model discrimination study, the initial value of the E-criterion is 6.41 x 10-6. Performing the dynamic optimisation, the optimal values for the design vector are presented diagrammatically in Figs. 6 and 7, giving a final value of the E-optimal criterion of 14.3 - - six orders of magnitude improvement in the predicted information content. The model predictions from the optimal design are shown in Fig. 6, where vertical arrows along the abscissa show the sampling points of the two response variables. Fig. 7(a) and (b) show the optimal input profiles for the dilution factor (Ul(t)) and feed substrate concentration (u2(t)), respectively. By simulating additional 'data' at these experimental conditions (variance as in §5.3), subsequent parameter estimation gives the results in Table 2. 0.3
40.0,
~
0.2
i:Z°.,
20.0
--1
L
,,(2,-.../,
\
~ ~ ~0.0. = =25.0. =.~ 20.0.
L,
;14•o1
.0 ¸
'; 15.0'
11.1 g 0.0
5
, ~, , ; . , . , , I0 IS 20 2S 30 3S 40 4S
5.0. 0.0 S
10
15
2O
25
3O
36
4O
rm.lkl
Fig. 7. The optimal input profiles for improving parameter precision.
2!°1
0•0 ~
0•O
' 5.0
',~ ~ ~I ,~ 10.0
15.0
20.0 25.0 Time [hi
30•0
, 35.0
~i 40.0
45.0
Fig. 6. The model predictions for improving parameter precision.
Fig. 5(a) and (b) show the optimal input profiles for the dilution factor (Ul(t)) and feed substrate concentration (u2(t)), respectively•
Table 2 Parameter estimates from the sequentially designed experiment for improving parameter precision Parameter
Estimate
Confidence interval
0, 02 03 04
0.311 0.254 0.523 0.046
+0.005 4-0.038 4-0.015 + 0.004
S.P. Asprey, S. Macchietto / Computers and Chemical Engineering 24 (2000) 1261-1267
As can be seen from these results, the confidence intervals have been reduced on average by 23.0% by using the designed experiments.
6. Concluding remarks In this paper, we have presented a general, systematic procedure to support the development and statistical verification of dynamic process models. Within this procedure, methods were presented to address several key aspects, such as structural identifiability and distinguishability testing, dynamic sensitivity analysis, and optimal design of dynamic experiments for both model discrimination and improving parameter precision. Within this framework, methods for accounting for parameter uncertainty were presented to provide experiment design robustness. Application of several of the tools was carried out using a typical biochemical application, demonstrating their effectiveness in aiding to build a dynamic process model.
1267
References Asprey, S. P., Macchietto, S., & Pantelides, C. C. (2000). Robust optimal designs for dynamic experiments. ADCHEM 2000. Pisa, Italy. Box, G. E. P., & Lucas, H. L. (1959). Design of experiments in non-linear situations. Biometrika, 46, 77-90. Espie, D. M., & Macchietto, S. (1989). The optimal design of dynamic experiments. American Institute of Chemical Engineering Journal, 35, 223-229. Hunter, W. G., & Reiner, A. M. (1965). Designs for discriminating between two rival models. Technometrics, 7, 307-323. Krrkel, S., Bauer, I., Bock, H. G., & Schlrder, J. P. (1999). A sequential approach for nonlinear optimum experimental design in DAE systems. Proceedings of the International Workshop on Scientific Computation in Chemical Engineering, vol. II. Pantelides, C. C. (2000). Personal communication. Stewart, W, E., Shon, Y., & Box, G. E. P. (1998). Discrimination and goodness of fit of multiresponse mechanistic models. American Institute of Chemical Engineering Journal, 44, 1404-1412. Vassiliadis, V. S., Sargent, R. W. H., & Pantelides, C. C. (1994). Solution of a class of multistage dynamic optimisation problems 1: problems without path constraints. Industrial & Engineering Chemistry Research, 33, 2111-2122. Walter, E. (1987). Identifiability of parametric models. Oxford: Pergamon Press. ZuUo, L. (1991). Computer aided design of experiments. An engineering approach. Ph.D. Thesis. University of London.
Computers and Chemical Engineering 24 (2000) 1261-1267
Computers &Chemical Engineering www.elsevier.com/locate/compchemeng
Statistical tools for optimal dynamic model building S.P. Asprey *, S. Macchietto Centre for Process Systems Engineering, Imperial College of Science, Technology and Medicine, Prince Consort Road, Londoln SW 7 2BY, UK
Abstract
A general, systematic procedure is presented to support the development and statistical verification of dynamic process models. Within this procedure, methods are presented to address several key aspects, such as structural identifiability and distinguishability testing, as well as optimal design of dynamic experiments for both model discrimination and improving parameter precision. A novel optimisation-based approach is introduced for testing of model structural identifiability and distinguishability, involving semi-infinite programming and max-rain problems. The design of dynamic experiments is cast as an optimal control problem within a framework that enables the calculation of optimal sampling points, experiment duration, fixed and variable external control profiles, and initial conditions of a dynamic experiment subject to general constraints on inputs and outputs. Within this framework, methods are presented to provide experiment design robustness, accounting for parameter uncertainty. The procedure is demonstrated through a practical biotechnology example. © 2000 Elsevier Science Ltd. All rights reserved. Keywords: Optimal experiment design; Dynamic modelling; Parameter estimation
1. Introduction
Detailed mathematical models are increasingly used by engineers in academia and industry to gain competitive advantage through such applications as modelbased process design, operations, and control. Thus, building and verifying high quality steady-state or dynamic, single or multiresponse mechanistic models of processing systems are key activities in process engineering. In this work, we present a general, systematic procedure to support the development and statistical verification of dynamic process models. The procedure comprises components that address several key aspects, such as structural identifiability and distinguishability testing, and optimal design of dynamic experiments for model discrimination and improving parameter precision. Experiment design is concerned with the following question: How does one adjust the time-varying controis, initial conditions or length of experiment for the experiment to generate the maximum amount of infor-
* Corresponding author. Tel.: + 44-20-75946604; fax: + 44-2075946606. E-mail address: [email protected] (S.P. Asprey)
mation for the purpose of determining the correct model from a set of candidates, and in turn, for estimating the parameters with greatest precision? For this, the problem is cast as an optimal control problem. Within this framework, one can calculate optimal sampiing points, experiment duration, fixed and variable external controls with input constraints, and initial conditions. To mathematically represent time-varying external controls to the process, we use the control vector parameterisation (CVP) technique, with either piecewise constant, or higher-order functions with 0th or higher-order continuity (Vassiliadis, Sargent & Pantelides, 1994). Here, we first present the overall procedure, followed by a brief discussion of current and new methods for each key component. Potential applications of these techniques in all engineering disciplines are abundant, including applications in chemical kinetics and reaction mechanism elucidation, polymer reaction engineering, physical properties estimation, and biochemical engineering. Throughout the paper, an engineering example using a typical semicontinuous bioreactor application is used to demonstrate the overall procedure and the power of the proposed techniques in their ability to reduce the quantity of experimental work required, while increasing the quality of the results.
0098-1354/00/$ - see front matter © 2000 Elsevier Science Ltd.'All rights reserved. PII: S0098-1354(00)00328-8
S.P. Asprey, S. Macchietto / Computers and Chemical Engineering 24 (2000) 1261-1267
1262
2. The modelling process When building models whose use is to explain an observed phenomenon, one uses a priori knowledge such as physical, chemical or biological 'laws' to propose (conceivably several) possible models. In each case, these laws dictate the model structure, and we may wish to know whether one or more such structures are adequate for the problem at hand. These models contain parameters that may have physical meaning, and we may wish to know if it is at all possible to determine their values, and, if so, to do so with maximum precision. In what follows, we therefore consider general deterministic models in the form of a set of (possibly mixed) differential and algebraic equations: f(~(t), x(t), u(t), w, 0, t) = 0 y(t) = g(x(t))
(1)
where x(t) is an ns-dimensional vector of time dependent state variables, u(t) is an no-dimensional vector of time-varying controls or inputs to the process, w is an n,.-dimensional vector of constant controls, 0 is a P-dimensional vector of model parameters to be determined within a continuous, realisable set 19, and y(t) is an M-dimensional vector of measured response variables that are a function of the state variables, x(t). We first present an overall model building strategy, depicted in Fig. 1. The first step to this strategy involves specifying one or more models to describe the process at hand. In Stage I, the strategy then uses these models to perform preliminary structural identifiability and distinguishability tests before any data are collected. These tests are used to determine whether or not the parameters in the mathematical formulation of the models can be uniquely identified, and whether or not the model structures can be distinguished from one another. In Stage II, the strategy involves designing experiments for model discrimination, followed by ..... .r..,.p..~.~.ry...A..,r.].y:!.~ ........ ~..~.fl .~ .,...a~....,r.n..,r. ....... .r.,.~..~.~..r...,~...L,., " ....
if* ©
E
parameter estimation and model adequacy checking to determine whether or not any of the proposed models can be rejected, with the final purpose of selecting the 'best' model from the given set. Once we are left with the 'best' model, we then go on to Stage III to design experiments for improving the precision of the parameters within that model, eventually arriving at a final, statistically-verified model formulation.
3. Preliminary analysis 3.1. Structural identifiability and distinguishability During the preliminary model analysis, we simply wish to know whether, given a candidate set of models and a proposed set of variables to be measured, we can in principle design an experiment such that the parameters of these models can be uniquely identified (identifiability), and whether the different model structures can be distinguished from one another (distinguishability). Past approaches for testing identifiability in nonlinear-in-input dynamic models find their roots in differential geometry, where, in particular, the model is approximated locally by its input-output mapping using functional expansions (Walter, 1987). Coefficients of these functional expansions are then evaluated at parameter values 0 and 0", and used to form systems of nonlinear algebraic equations, whose unique solution 0 = 0 " ensures identifiability. This approach has its advantages and disadvantages; its advantages being (i) it is a definitive test (answers the test with yes/no), and (ii) the methods are not affected by scaling problems. Disadvantages of the methods are (i) they are limited to models of special forms (e.g. nonlinear control-affine models), and (ii) they are limited to models with dimensions (ns+ P)_< 10, very small for modelling realistic process systems. To avoid these disadvantages (in particular disadvantage (ii)), we propose a novel optimisation-based approach, defined in the following (Pantelides, 2000).
Definition: Identifiability. A model giving output trajectory y(0, u(t), t) is globally identifiable if, for any two parameter sets, 0 c O and 0"~O*, a time horizon of interest, t~[0, tA, for all system inputs u(t)eU, and the same initial conditions y(t = 0), the global maximum: Test II:
~z=
max
0~0,0"~0"
(0-0*)zW0(0-0*)
(2a)
subject to: (y(u(t), O) -- y(u(t), O*))rWy(y(u(t), O) -- y(u(t), 0")) Fig. 1. The overall model-building scheme.
x dt < eyV u(t)eU
S.P. Asprey, S. Macchietto / Computers and Chemical Engineering 24 (2000) 1261-1267 12' 1' 0.8,
0.2, I ti=100
I
0 9a
40
60
80
I00
I~
Fig. 2. Structural identifiabilityand distinguishabilityillustration.
f(~(t), y(t), x(t), u(t), w, 0, t) = 0 (the model equations)
0~ < 0 , < 0~v 07L < 0 * < 0 T v
i=1 ..... P (2b)
Definition: Distinguishability. Two models giving output trajectories y(O, u(t), t) and y*(0*, u(t), t), respectively, are distinguishable if for any 0~O, 0*eO*, re[O, t/] for all u(t)eU, and the same initial conditions y(t = 0), the global maximum: Test DI:
fT min
u(t)EU0~0,0"~0" 0
(yo)(t) - y(2)(t)) r
x W(yo)(t ) - y{2)(t)) dt
(3a)
fo)(£(t), y(t), x(t), u(t), w, 0, t) = 0 f(2)(~,(t), y(t), x(t), u(t), w, 0", t) = 0
i=l,...,P
0 * L < 0 * < 0 *t~ i = l . . . . . P
4.1. Model d&crimination The problem of model discrimination occurs when ArM rival models have been proposed to describe a system and it is not certain which of the ArM models is 'best'. For any experiment, it is expected that the 'best' model will provide the most accurate prediction of the observed measurements. Experimental designs optimised for model discrimination are termed T-optimal, and are based on maximising an overall measure of the divergence between the model predictions. Very little past work has been presented for discriminating between rival dynamic models. Espie and Macchietto (1989) were the first to attempt an extension of the Hunter and Reiner (1965) criterion to the dynamic case, where the summation over the data points in the original formulation was replaced by an integral over a pre-specified time interval [to, t/]: ft/~ Um ^^ q~r = arg max ~ Tt,f(g~,O, 0', t) dt (4) r,,,,(~,
subject to:
Of
responses of the two models, over the time horizon of interest, If such a value is sufficiently large, the two models may be deemed distinguishable. Again, returning to our illustration of Fig. 2, we wish to ensure that when the same input u(t) is applied to two different models with parameter vectors 0 and 0", respectively, and the same initial conditions, the two predicted model output trajectories do not overlie. This optimisation approach is not limited to any particular model form or size; however, we are now faced with solving semi-infinite or m a x - m i n dynamic optimisation problems. We must also choose a priori tolerances for each test (E% and ~ , best chosen to be commensurate with the integral o~ the expected variance of experimental measurements). A more complete derivation of the optimisation-based approach to testing identifiability and distinguishability, as well as solution methods will be given in a future paper.
4. Design of dynamic experiments
i=l ..... P
gives @~< %. Here e%and ey are arbitrarily small positive values and W0 and Wy are weighting matrices. In words, @z is the largest distance between two parameter vectors that still give (essentially) the same response trajectories. If such a distance is arbitrarily small, the model can be deemed globally identifiable. In our illustration of Fig. 2, we wish to ensure that when the same input u(t) is applied to the same model, at the same initial conditions, but with different parameter vectors 0 and 0* (such that 0 ~ 0"), the two predicted model output trajectories do not overlie.
@D = max
1263
(3b)
gives @D > ~ o for some ~ o > 0. In words, @D is the largest value of the minimum distance between the
~, ~', t)
=
(~,(e, ~, t)
- ~,,(to, g', t)) ~
× W(~t(q), 0, t) -- ~r(q), 0', t))
(5)
and embedded in an optimal control framework. Here, is the vector of the best available estimates of the model parameters, and (p is the vector of experiment design variables. However, by using a continuous function (Eq. (4)), the meaning of the sampling time of the responses is lost. In our newly proposed formulation, we keep the concept of sampling time by using the Hunter-Reiner formulation as a starting point for the design of dynamic experiments, where now we use the notion of discrete sampling points, tk, k = 1 , . . . , nsp of each of the response measurements:
S.P. Asprey, S. Macchietto / Computers and Chemical Engineering 24 (2000) 1261-1267
1264
arg max
~1[~T =
NM
NM
~,
~
T4r(~O,O, 0', tk)
(6)
~p~O k = l l = l 1'=l--1
We embed this criterion within a dynamic optimisation framework, described in detail in Asprey, Macchietto and Pantelides (2000). The design vector, q~, now comprises the discrete sampling points as well as any time-varying and time-invariant inputs, and variable initial conditions. To account for uncertainty in the point estimates of the parameters, and thus variance in the predictions of the models, new extensions of (6) could be formulated based on the Bayesian posterior density approach used for nonlinear steady-state process models (Stewart, Shon & Box, 1998), or by formulating either an expected value criterion: qJTEV
= arg max q~
f ,~= NM
E
"~ ~
O~O,O'EO' k
NM
1 l=1 l'=1+1
--}
min ~[V(0, ¢p, ws)]. Cp,W
Tu,(~p, 0, 0', tk)
(7) or a worst-case criterion:
tPrwc = arg max ~
min
NM
NM
Y',
~
However, this formulation is restricted to a diagonal variance-covariance matrix of the parameters, an assumption that is almost always false. Again, the continuous function (Eq. (11)) approach does not explicitly allow one to optimise discrete sampling times of the system responses. More recent work has been carried out by KSrkel, Bauer, Bock and Schl6der (1999) on sequential experiment design for parameter precision improvement in models described by DAEs. They take a similar optimal control approach, using piecewise constant or piecewise linear functions for parameterising the time-varying controls, but instead use a metric (@[] = determinant, D-optimal; maximum eigenvalue, E-optimal, etc.) of the parameter variance-covariance matrix as an objective function:
Tt,r(~p, 0, @, tk)
OeO,O'~O' k = 1 1=1 l ' = l + 1
(8)
For selecting optimal sampling times of the system responses, they use binary decision variables, ws, as weights to select from a set of equally spaced points along the time horizon of interest. Much of our work to date is based on a unique criterion put forth by Zullo (1991), who defined a discrete version of the dynamic information matrix for experiment design for the improvement of parameter precision, defined by:
4.2. Improving parameter precision
M
=
The aim of experiment design for improved parameter precision is to decrease the size of the inference regions of each of the parameters in any given model. This is equivalent to making the elements of the parameter variance-covariance matrix 'small'. Thus, we use the marginal posterior density covariance, given by (Box & Lucas, 1959): V(O, q~) :'~
[~ f r
Vl'~'rs~T~', ~'~
]--1
=
[ Srl
Sr2
Qr = ISrl
(10)
~r{~ is the rsth element of the inverse of the estimate of the variance-covariance matrix ~y = cov(y,, y~). Espie and Macchietto (1989) were the first to attempt an extension of these concepts to the dynamic case for systems described by DAEs, formulating the design of experiments as an optimal control problem and suggesting a criterion based on maximising the weighted sum of the squared dynamic sensitivity coefficients: ~ wi/(S,y)2 dt
qJ~ ¢l~ .it 0 j = 1 i ~ 1
""
Srp
] ~
S,2 "'" S~e[ ~
1 S,2 •
80
q~ = arg max
where tp is the vector of experiment design variables, as above. The matrix Qr is the matrix of dynamic sensitivity coefficients of the rth response variable in the model computed at each of the n~e sampling points:
ls=1
8y,(cp,0, t) o=~
(12)
r=ls=l
sampling point sampling point
t~p1 t~p2
sampling point
tsp,sp
(9)
where Sr is the vector of partial derivatives of the rth response variable in the model with respect to the parameters 0 calculated at all experimental points: S,
M
E
(11)
"'" ,
Sr e J
~
(13) Extensions of this criterion to account for the uncertainty in the point estimates of the parameters are made in Asprey et al. (2000), in which the D-optimal criterion is embedded within a robust optimisation framework, using either an expected value approach, based upon a multivariate normal probability distribution function (pdf) describing the prior uncertainty in the parameters, 0 ,,~ N(O, 2~):
q~ev = arg max E {[Mz(0, q~)]} q~
(14)
0~
where: oE { I M / ( o,
,p)I}=
IM/(0, ¢p)
S.P. Asprey, S. Macchietto / Computers and Chemical Engineering 24 (2000) 1261-1267
or a worst-case approach over an admissible domain O: • R = arg max min {[MI(0, tp)[} tpe~,
(15)
5. A process engineering example
Here, an example is presented to illustrate the application of the dynamic model building tools presented in Sections 2-4. Consider a fed-batch reactor (see Fig. 3) in which the fermentation of baker's yeast is carried out. Within this process, there are two time-varying controls (ul, Uz) and two measured concentrations
1265
5.1. Identifiability o f Mz(O, u, t)
Global structural identifiability of MI(O, u, t) can be tested via Test I1 and using solution techniques for such semi-infinite programming problems (see Asprey et al., 2000). The solution, using 0~e[0.01, 0.98] Vi, a five-element piece-wise constant CVP for the process inputs and a time horizon of interest of 40 h, was found to be ~ y = 1.166 x 10 -6 ( < g.a,t --- 1.0 x 10-5; gy :- 4.0 x 10-3), n o t significantly different from zero and thus indicating the global identifiability of 34/(0, u, t). 5.2. Distinguishability o f MI(O, u, t) and MIr(O, u, t)
( X l , X2).
We propose two models to describe the system. Two differential equations, common to both models, describe the component balances within the reactor: dXl dt = (r - Ua - 04)xl dx--3 = - rx--2+ ul(u2 - x2)
dt
(16)
03
The first model assumes Monod-type kinetics for substrate consumption, giving a third model equation: MI(0, u, t):
01X2
r=
02 + x2
(17)
The second model assumes Canoid-type kinetics and thus, instead of Eq. (17), gives the model form: MII(O, u, t):
01x2
r=
(18)
0 2 X 1 -4- X 2
The experimental conditions that characterise a particular experiment are: (i) the initial biomass concentration (or inoculation), x °, with range 1-10 g 1-~; (ii) the dilution factor, ua, with range 0.05-0.20 h-1; (iii) the substrate concentration in the feed, u2, with range 5-35 g 1-1. The initial substrate concentration, x2(0), is always 0.01 g 1-~ and cannot be manipulated for experiment design purposes. Both x~ and x2 can be measured during the experiment (y = [xb x2]r).
nip~uts
a~-v~n¢"'~,_
Global structural distinguishability of Mz(0, u, t) and Mn(O, u, t) can be tested by reformulating the m a x -
min Test D1 problem to a semi-infinite programming problem and using the techniques mentioned above. The resulting solution of the distinguishability test using the same conditions and CVP as above was found to be ~ D = 3.392 x 10 -1, giving a maximum larger than e~o.= 4.0 × 10 - 3 and thus confirming the distinguishabdity between MI(0, u, t) and MzI(O, u, t). 5.3. Designing a dynamic experiment for model discrimination
Using models 3/1i(0, u, t) and MI(0, u, t), we begin by assuming availability of initial parameter estimates 01 = [0.30, 0.25, 0.56, 0.02] r and On = [0.30, 0.03, 0.55, 0.03] r, respectively. A dynamic experiment design was then carried out for model discrimination using Eq. (6) and piece-wise constant control inputs, with a maximum experiment duration of 40 h. An a priori selection of the number of sampling points, along with the number of switching times/constant levels of the input profiles is required. We arbitrarily chose ten sampling times (nsp = 10) and five time intervals for the control vector parameterisation (Gwi = 5, i = l, " ' " ' n, where, for the • bloreactor models nu 2). We start the optimisation with an experiment lasting I0 h, with hourly sampling points and constant control inputs at 0.2 and 15.0), giving an initial divergence of 1.98. We impose the constraints (cf. Asprey et al., 2000): 1.0
_
5.0
Fig. 3. A fed-batch reactor.
i =
1. . . . ,nsv
i = 1 . . . . . G;
j = 1. . . . . Gw.~ The optimal values for the design vector maximising criterion (6) are presented diagrammatically in Figs. 4 and 5, giving an optimal divergence of 1.19 x 105. The model predictions from the optimal design are shown in Fig. 4. The optimal experiment lasts 40 h and vertical arrows along the abscissa show the optimal sampling points of the two response variables (both measured at the same sampling points).
S.P. Asprey, S. Macchietto / Computers and Chemical Engineering 24 (2000) 1261-1267
1266
By simulating additional 'data' at these experimental conditions, with a homoscedastic normal error variance of 0.04 for both responses, subsequent parameter estimation gives the results in Table 1.By inspection of the multiresponse Z 2 lack-of-fit analyses shown in Table 1, we can conclude that model MI(0, u, t) is better suited to describe the observed experimental data, although, statistically, both models describe the observed data very well.
111.0
~
10•0
'
14,0
'
12•0
'
0.0
/'0.0
YOI(21- . , . . ~ . 5.0
,.2.0
.
. 10.0
.
I~
tt
. 15.0
20.0
25.0
Time
30.0
35.0
45.0
40.0
[11]
5.4. Designing a dynamic experiment for improving parameter precision
Fig. 4. The model predictions for model discrimination.' 0.3.
36.0,
£0.2-
g 3S.0,
LU
li.o 0.1" ! 0 . 1 - ~ 0.0 10.0
0.0
~ 31.0'
,
_, : _,
20.0
30.0
40.0
50.0
300 0.0
10.0
TmN
- - 20.0 30.0 llmet~l
40.0
50,0
Fig. 5. The optimal input profiles for model discrimination.
Table 1 Parameter estimates for MI(0, u, t) and Mn(0, u, t) after one designed experiment Parameter
Estimate
Confidence interval
0= 02 03
Mz 0.325 0.267 0.559 0.055
Mu 0.303 0.011 0.549 0.049
Mz 40.012 +0.075 __40.018 + 0.006
0.58 1.96
Z2ef= 26.3 Z2ef= 26.3
04
LOF (M/) LOF (Mzz)
MI! _+0.011 +0.009 _+0.023 _+0.011
We next turn to improving the precision in the parameter estimates in MI(O, u, t). Using the estimates shown in Table 1, an experiment design for improving parameter precision using the E-optimality criterion with Eq. (12) was performed. Using the same initial values, constraints, and representation of the time-varying inputs as in the model discrimination study, the initial value of the E-criterion is 6.41 x 10-6. Performing the dynamic optimisation, the optimal values for the design vector are presented diagrammatically in Figs. 6 and 7, giving a final value of the E-optimal criterion of 14.3 - - six orders of magnitude improvement in the predicted information content. The model predictions from the optimal design are shown in Fig. 6, where vertical arrows along the abscissa show the sampling points of the two response variables. Fig. 7(a) and (b) show the optimal input profiles for the dilution factor (Ul(t)) and feed substrate concentration (u2(t)), respectively. By simulating additional 'data' at these experimental conditions (variance as in §5.3), subsequent parameter estimation gives the results in Table 2. 0.3
40.0,
~
0.2
i:Z°.,
20.0
--1
L
,,(2,-.../,
\
~ ~ ~0.0. = =25.0. =.~ 20.0.
L,
;14•o1
.0 ¸
'; 15.0'
11.1 g 0.0
5
, ~, , ; . , . , , I0 IS 20 2S 30 3S 40 4S
5.0. 0.0 S
10
15
2O
25
3O
36
4O
rm.lkl
Fig. 7. The optimal input profiles for improving parameter precision.
2!°1
0•0 ~
0•O
' 5.0
',~ ~ ~I ,~ 10.0
15.0
20.0 25.0 Time [hi
30•0
, 35.0
~i 40.0
45.0
Fig. 6. The model predictions for improving parameter precision.
Fig. 5(a) and (b) show the optimal input profiles for the dilution factor (Ul(t)) and feed substrate concentration (u2(t)), respectively•
Table 2 Parameter estimates from the sequentially designed experiment for improving parameter precision Parameter
Estimate
Confidence interval
0, 02 03 04
0.311 0.254 0.523 0.046
+0.005 4-0.038 4-0.015 + 0.004
S.P. Asprey, S. Macchietto / Computers and Chemical Engineering 24 (2000) 1261-1267
As can be seen from these results, the confidence intervals have been reduced on average by 23.0% by using the designed experiments.
6. Concluding remarks In this paper, we have presented a general, systematic procedure to support the development and statistical verification of dynamic process models. Within this procedure, methods were presented to address several key aspects, such as structural identifiability and distinguishability testing, dynamic sensitivity analysis, and optimal design of dynamic experiments for both model discrimination and improving parameter precision. Within this framework, methods for accounting for parameter uncertainty were presented to provide experiment design robustness. Application of several of the tools was carried out using a typical biochemical application, demonstrating their effectiveness in aiding to build a dynamic process model.
1267
References Asprey, S. P., Macchietto, S., & Pantelides, C. C. (2000). Robust optimal designs for dynamic experiments. ADCHEM 2000. Pisa, Italy. Box, G. E. P., & Lucas, H. L. (1959). Design of experiments in non-linear situations. Biometrika, 46, 77-90. Espie, D. M., & Macchietto, S. (1989). The optimal design of dynamic experiments. American Institute of Chemical Engineering Journal, 35, 223-229. Hunter, W. G., & Reiner, A. M. (1965). Designs for discriminating between two rival models. Technometrics, 7, 307-323. Krrkel, S., Bauer, I., Bock, H. G., & Schlrder, J. P. (1999). A sequential approach for nonlinear optimum experimental design in DAE systems. Proceedings of the International Workshop on Scientific Computation in Chemical Engineering, vol. II. Pantelides, C. C. (2000). Personal communication. Stewart, W, E., Shon, Y., & Box, G. E. P. (1998). Discrimination and goodness of fit of multiresponse mechanistic models. American Institute of Chemical Engineering Journal, 44, 1404-1412. Vassiliadis, V. S., Sargent, R. W. H., & Pantelides, C. C. (1994). Solution of a class of multistage dynamic optimisation problems 1: problems without path constraints. Industrial & Engineering Chemistry Research, 33, 2111-2122. Walter, E. (1987). Identifiability of parametric models. Oxford: Pergamon Press. ZuUo, L. (1991). Computer aided design of experiments. An engineering approach. Ph.D. Thesis. University of London.