Невозможно использовать функции из пакетов R, fda и pCODE вместе из-за ошибки: "аргумент" раз "отсутствует, без значения по умолчанию" - PullRequest
Новая
       146

Невозможно использовать функции из пакетов R, fda и pCODE вместе из-за ошибки: "аргумент" раз "отсутствует, без значения по умолчанию"

0 голосов
/ 06 августа 2020

В настоящее время я изучаю уравнения Tank Reactor из Ramsay et al и wi sh, чтобы найти матрицу дисперсии, которую он использовал. Я использую пакет R, fda, ​​в частности, функции CSTR2 и CSTR2in в сочетании с пакетом pCODE, используя функцию deltavar, чтобы найти матрицу дисперсии. Я установил ODE аналогично тому, как Рамзи установил его в своем пакете fda: 

library(pCODE)
library(deSolve)
library(fda)
library(MASS)
library(pracma)
library(Hmisc)

#Here fitstruct is our parameter list
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
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
times = seq(0, Tlim, delta)#

#These are the 5 inputs we are interested in
CSTR_Inputs <- CSTR2in(times, 'all.cool.step')

#  set constants for ODE solver

#  cool condition solution
#  initial conditions

C_initial = 1.5965#  initial concentration in kmol per cubic meter
T_initial = 341.3754# initial temperature in deg K
initial_values = c(Conc = C_initial, Temp=T_initial)

#  load cool input into fitstruct
fitstruct[["Tcin"]] = CSTR_Inputs[, "Tcin"];

CSTR <- function (times, state, parameters) 
{
  tau = 1
  CSTR.in <- CSTR2in(times, parameters$condition, tau)
  fitstruct <- parameters$fitstruct
  kref = fitstruct$kref
  EoverR = fitstruct$EoverR
  a = fitstruct$a
  b = fitstruct$b
  V = fitstruct$V
  FoverV = CSTR.in[, "F."]/V
  aFc2bp = a * CSTR.in[, "Fc"]^b
  const1 = V * fitstruct$rho * fitstruct$Cp
  const2 = 2 * fitstruct$rhoc * fitstruct$Cpc
  betaTTcin = (CSTR.in[, "Fc"] * aFc2bp/(const1 * (CSTR.in[, 
                                                           "Fc"] + aFc2bp/const2)))
  betaTT0 = FoverV
  {
    if (length(state) == 2) {
      Ci = state[1]
      Ti = state[2]
    }
    else {
      Ci = state[, 1]
      Ti = state[, 2]
    }
  }
  Tdif = 1/as.vector(Ti) - 1/fitstruct$Tref
  betaCC = kref * exp(-10000 * EoverR * Tdif)
  betaTC = (-fitstruct$delH/const1) * betaCC
  betaTT = FoverV + betaTTcin
  DC = (-(FoverV + betaCC) * as.vector(Ci) + (CSTR.in[, "F."]/V) * 
          CSTR.in[, "CA0"])
  DT = (-betaTT * Ti + betaTC * Ci + (CSTR.in[, "F."]/V) * 
          CSTR.in[, "T0"] + betaTTcin * CSTR.in[, "Tcin"])
  x <- cbind(Conc = DC, Temp = DT)
  list(x)
}

И использовал pCODE следующим образом: 

model.par<-list(fitstruct=fitstruct, condition='all.cool.step')
#fda recommends using lsoda instead of ode function
CSTR_Solution <- ode(y=initial_values, times=times, func=CSTR,
                        parms=model.par)

n_obs <- length(times)

#setting up observed data
observations <- matrix(NA, nrow = length(times),ncol =3)
observations[,1] <- times
observations[,2] <- CSTR_Solution[,2] + rnorm(n = n_obs, mean = 0 , sd = 0.0223)
observations[,3] <- CSTR_Solution[,3] + rnorm(n = n_obs, mean = 0 , sd = 0.79)

# 100 knots equally spaced within [0,20] as 401 knots are too much in this case.
# Nested optimization usually takes too long.
knots <- seq(from=0, to=Tlim, length.out = 101)

#order of basis functions
norder <- 3

#number of basis funtions
nbasis <- length(knots) + norder - 2

#creating Bspline basis
basis  <- create.bspline.basis(c(0,Tlim),nbasis,norder,breaks = knots)

deltavar(data = observations[,2:3], time = times, ode.model = CSTR,
         par.initial = c(0.8,0.4,1.7), par.names = c('EoverR','kref','a'),state.names = c('C','T'),
         basis.list = list(basis,basis), lambda = 1e2,
         stepsize = 1e-3,y_stepsize = 1e-3)

Однако, появляется следующая ошибка: 

Error in CSTR2in(times, parameters$condition, tau) : 
  argument "times" is missing, with no default

Я не совсем уверен, как это исправить, вот трассировка, если это поможет. 

14: CSTR2in(times, parameters$condition, tau) at #4
13: derive.model(state = x, parameters = temp)
12: unlist(derive.model(state = x, parameters = temp))
11: FUN(newX[, i], ...)
10: apply(Xt, 1, temp_fun)
9: t(apply(Xt, 1, temp_fun))
8: fct(x, ...)
7: fun(x0)
6: nls_optimize.inner(innerobj_multi, basis.initial, ode.par = ode.parameter, 
       derive.model = derivative.model, options = list(maxeval = 50, 
           tolx = 1e-06, tolg = 1e-06), input = inner.input)
5: fct(x, ...)
4: fun(x0)
3: nls_optimize(options = list(maxeval = con.now$maxeval, tau = con.now$tau), 
       outterobj_multi_nls, par.initial, basis.initial = unlist(initial_coef), 
       derivative.model = ode.model, inner.input = inner.input, 
       verbal = 1)
2: pcode(data, time, ode.model, par.names, state.names, likelihood.fun, 
       par.initial, basis.list, lambda = lambda, controls = con.now)
1: deltavar(data = observations[, 2:3], time = times, ode.model = CSTR, 
       par.initial = c(0.8, 0.4, 1.7), par.names = c("EoverR", "kref", 
           "a"), state.names = c("C", "T"), basis.list = list(basis, 
           basis), lambda = 100, stepsize = 0.001, y_stepsize = 0.001)
...