Data Engineering and BI courses are free this week!
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)
CSTRres0(kref=NULL, EoverR=NULL, a=NULL, b=NULL, gradwrd=FALSE)
CSTRsse(par, datstruct, fitstruct, CSTRbasis, lambda)
Sres = a matrix giving the residuals between observed and predicted datstruct[["y"]] divided by sqrt(datstruct[[c("Cwt", "Twt")]]) so the result is dimensionless. dim(Sres) = dim(datstruct[["y"]]). Thus, if datstruct[["y"]] has only one column, 'Sres' has only one column.
Lres = a matrix with two columns giving the difference between
left and right hand sides of the CSTR differential equation at
all the quadrature points. dim(Lres) = c(nquad, 2).
}
DSres = a matrix with one row for each element of res[["Sres"]] and two columns for each basis function.
DLres = a matrix with two rows for each quadrature point and two columns for each basis function.
If gradwrd=FALSE, this component is not present. }
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)
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 two functions CSTRres and CSTRres0 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 'odesolve' 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 Gauss-Newton optimization to estimates the coefficients of CSTRbasis to minimize the weighted sum of squares of residuals returned by CSTRfitLS.
CSTRres and CSTRres0 provide alternative interfaces between 'nls' and 'CSTRfn'. Both get the parameters to be estimated via their 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'. (The version of 'nls' in the 'stats' package in R 2.5.0 required 'data' to be a data.frame, not merely a list. The version of 'nls' in 'fda' does not require this.) By contrast, CSTRres0 uses 'get' to obtain these other arguments as .datstruct, .fitstruct, .CSTRbasis and .lambda. Thus, before calling nls(~CSTRres0(...), ...), the values of these arguments must be assigned to .datstruct, .fitstruct, .CSTRbasis and .lambda. If 'nls' is called from the global environment, it will look for these objects in the global environment.
CSTRres0 has one feature absent from CSTRres: If a variable .CSTRres0.trace is available, it is assumed to be a matrix with columns kref, EoverR, a, b, SSE, plus all residuals. These numbers are rbinded as an additional row of this matrix. This is provided to help diagnose a problem were 'nls' was terminating with "step factor ... reduced below 'minFactor'", facilitating the comparison between R and Matlab for the precise sets of parameter values tested by 'nls'.
CSTRsse computes sum of squares of residuals for use with optim or nlminb.
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. 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 899-902.
##
## 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(odesolve)) {
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.Run the code above in your browser using DataLab