Functions for solving the Continuously Stirred Tank Reactor (CSTR) Ordinary Differential Equations (ODEs). A solution for observations where metrology error is assumed to be negligible can be obtained via lsoda(y, Time, CSTR2, parms); CSTR2 calls CSTR2in. When metrology error can not be ignored, use CSTRfn (which calls CSTRfitLS). To estimate parameters in the CSTR differential equation system (kref, EoverR, a, and / or b), pass CSTRres to nls. If nls fails to converge, first use optim or nlminb with CSTRsse, then pass the estimates to nls.
1 2 3 4 5 6 7 8 9 10 11 12  CSTR2in(Time, condition =
c('all.cool.step', 'all.hot.step', 'all.hot.ramp', 'all.cool.ramp',
'Tc.hot.exponential', 'Tc.cool.exponential', 'Tc.hot.ramp',
'Tc.cool.ramp', 'Tc.hot.step', 'Tc.cool.step'),
tau=1)
CSTR2(Time, y, parms)
CSTRfitLS(coef, datstruct, fitstruct, lambda, gradwrd=FALSE)
CSTRfn(parvec, datstruct, fitstruct, CSTRbasis, lambda, gradwrd=TRUE)
CSTRres(kref=NULL, EoverR=NULL, a=NULL, b=NULL,
datstruct, fitstruct, CSTRbasis, lambda, gradwrd=FALSE)
CSTRsse(par, datstruct, fitstruct, CSTRbasis, lambda)

Time 
The time(s) for which computation(s) are desired 
condition 
a character string with the name of one of ten preprogrammed input scenarios. 
tau 
time for exponential decay of exp(1) under condition = 'Tc.hot.exponential' or 'Tc.cool.exponential'; ignored for other values of 'condition'. 
y 
Either a vector of length 2 or a matrix with 2 columns giving the observation(s) on Concentration and Temperature for which computation(s) are desired 
parms 
a list of CSTR model parameters passed via the lsoda 'parms' argument. This list consists of the following 3 components:

coef 
a matrix with one row for each basis function in fitstruct and columns c("Conc", "Temp") or a vector form of such a matrix. 
datstruct 
a list describing the structure of the data. CSTRfitLS uses the following components:

fitstruct 
a list with 14 components:

lambda 
a 2vector of rate parameters 'lambdaC' and 'lambdaT'. 
gradwrd 
a logical scalar TRUE if the gradient is to be returned as well as the residuals matrix. 
parvec, par 
initial values for the parameters specified by fitstruct[[ "estimate"]] to be estimated. 
CSTRbasis 
Quadrature basis returned by 'quadset'. 
kref, EoverR, a, b 
the kref, EoverR, a, and b coefficients of the CSTR model as individual arguments of CSTRres to support using 'nls' with the CSTR model. Those actually provided by name will be estimated; the others will be taken from '.fitstruct'; see details. 
Ramsay et al. (2007) considers the following differential equation system for a continuously stirred tank reactor (CSTR):
dC/dt = (betaCC(T, F.in)*C + F.in*C.in)
dT/dt = (betaTT(Fcvec, F.in)*T + betaTC(T, F.in)*C + alpha(Fcvec)*T.co)
where
betaCC(T, F.in) = kref*exp(1e4*EoverR*(1/T  1/Tref)) + F.in
betaTT(Fcvec, F.in) = alpha(Fcvec) + F.in
betaTC(T, F.in) = (delH/(rho*Cp))*betaCC(T, F.in)
alpha(Fcvec) = (a*Fcvec^(b+1) / (K1*(Fcvec + K2*Fcvec^b)))
K1 = V*rho*Cp
K2 = 1/(2*rhoc*Cpc)
The four functions CSTR2in, CSTR2, CSTRfitLS, and CSTRfn compute coefficients of basis vectors for two different solutions to this set of differential equations. Functions CSTR2in and CSTR2 work with 'lsoda' to provide a solution to this system of equations. Functions CSTSRitLS and CSTRfn are used to estimate parameters to fit this differential equation system to noisy data. These solutions are conditioned on specified values for kref, EoverR, a, and b. The other function, CSTRres, support estimation of these parameters using 'nls'.
CSTR2in translates a character string 'condition' into a data.frame containing system inputs for which the reaction of the system is desired. CSTR2 calls CSTR2in and then computes the corresponding predicted first derivatives of CSTR system outputs according to the right hand side of the system equations. CSTR2 can be called by 'lsoda' in the 'deSolve' package to actually solve the system of equations. To solve the CSTR equations for another set of inputs, the easiest modification might be to change CSTR2in to return the desired inputs. Another alternative would be to add an argument 'input.data.frame' that would be used in place of CSTR2in when present.
CSTRfitLS computes standardized residuals for systems outputs Conc, Temp or both as specified by fitstruct[["fit"]], a logical vector of length 2. The standardization is sqrt(datstruct[["Cwt"]]) and / or sqrt(datstruct[["Twt"]]) for Conc and Temp, respectively. CSTRfitLS also returns standardized deviations from the predicted first derivatives for Conc and Temp.
CSTRfn uses a GaussNewton optimization to estimates the coefficients of CSTRbasis to minimize the weighted sum of squares of residuals returned by CSTRfitLS.
CSTRres provides an interface between 'nls' and 'CSTRfn'. It gets the parameters to be estimated via the official function arguments, kref, EoverR, a, and / or b. The subset of these paramters to estimate must be specified both directly in the function call to 'nls' and indirectly via fitstruct[["estimate"]]. CSTRres gets the other CSTRfn arguments (datstruct, fitstruct, CSTRbasis, and lambda) via the 'data' argument of 'nls'.
CSTRsse computes sum of squares of residuals for use with optim or nlminb.
CSTR2in 
a matrix with number of rows = length(Time) and columns for F., CA0, T0, Tcin, and Fc. This gives the inputs to the CSTR simulation for the chosen 'condition'. 
CSTR2 
a list with one component being a matrix with number of rows = length(tobs) and 2 columns giving the first derivatives of Conc and Temp according to the right hand side of the differential equation. CSTR2 calls CSTR2in to get its inputs. 
CSTRfitLS 
a list with one or two components as follows:

CSTRfn 
a list with five components:

CSTRres 
the 'res' component of CSTRfd(...) as a column vector. This allows us to use 'nls' with the CSTR model. This can be especially useful as 'nls' has several helper functions to facilitate evaluating goodness of fit and and uncertainty in parameter estimates. 
CSTRsse 
sum(res*res) from CSTRfd(...). This allows us to use 'optim' or 'nlminb' with the CSTR model. This can also be used to obtain starting values for 'nls' in cases where 'nls' fails to converge from the initiall provided starting values. Apart from 'par', the other arguments 'datstruct', 'fitstruct', 'CSTRbasis', and 'lambda', must be passed via '...' in 'optim' or 'nlminb'. 
Ramsay, J. O., Hooker, G., Cao, J. and Campbell, D. (2007) Parameter estimation for differential equations: A generalized smoothing approach (with discussion). Journal of the Royal Statistical Society, Series B, 69, 741796.
Ramsay, J. O., and Silverman, B. W. (2006) Functional Data Analysis, 2nd ed., Springer, New York.
Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.
lsoda
nls
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66  ###
###
### 1. lsoda(y, times, func=CSTR2, parms=...)
###
###
# The system of two nonlinear equations has five forcing or
# input functions.
# These equations are taken from
# Marlin, T. E. (2000) Process Control, 2nd Edition, McGraw Hill,
# pages 899902.
##
## Set up the problem
##
fitstruct < list(V = 1.0,# volume in cubic meters
Cp = 1.0,# concentration in cal/(g.K)
rho = 1.0,# density in grams per cubic meter
delH = 130.0,# cal/kmol
Cpc = 1.0,# concentration in cal/(g.K)
rhoc = 1.0,# cal/kmol
Tref = 350)# reference temperature
# store true values of known parameters
EoverRtru = 0.83301# E/R in units K/1e4
kreftru = 0.4610 # reference value
atru = 1.678# a in units (cal/min)/K/1e6
btru = 0.5# dimensionless exponent
#% enter these parameter values into fitstruct
fitstruct[["kref"]] = kreftru#
fitstruct[["EoverR"]] = EoverRtru# kref = 0.4610
fitstruct[["a"]] = atru# a in units (cal/min)/K/1e6
fitstruct[["b"]] = btru# dimensionless exponent
Tlim = 64# reaction observed over interval [0, Tlim]
delta = 1/12# observe every five seconds
tspan = seq(0, Tlim, delta)#
coolStepInput < CSTR2in(tspan, 'all.cool.step')
# set constants for ODE solver
# cool condition solution
# initial conditions
Cinit.cool = 1.5965# initial concentration in kmol per cubic meter
Tinit.cool = 341.3754# initial temperature in deg K
yinit = c(Conc = Cinit.cool, Temp=Tinit.cool)
# load cool input into fitstruct
fitstruct[["Tcin"]] = coolStepInput[, "Tcin"];
# solve differential equation with true parameter values
if (require(deSolve)) {
coolStepSoln < lsoda(y=yinit, times=tspan, func=CSTR2,
parms=list(fitstruct=fitstruct, condition='all.cool.step', Tlim=Tlim) )
}
###
###
### 2. CSTRfn
###
###
# See the script in '~R\library\fda\scripts\CSTR\CSTR_demo.R'
# for more examples.

Questions? Problems? Suggestions? Tweet to @rdrrHQ or email at ian@mutexlabs.com.
All documentation is copyright its authors; we didn't write any of that.