如何使用rtop包中的sceua方法优化模型参数，R语言

Question

我是一名对水文学和径流模拟感兴趣的学生。我在 R 中使用 topmodel 包和 huagrahuma 数据集，我想通过在 rtop 包中使用 sceua 方法来优化它的参数。我写了一个代码，使用“rtop package | R Documentation”中的示例作为参考。但是，我在优化参数时出错。它说“错误在 if (sum(mapply(FUN = function(x, y, z) max(y - x, x - z, 0), : missing value where TRUE/FALSE needed”。起初，我怀疑错误是由Qobs具有的NA值触发的。所以我将NA转换为数值，但再次出现相同的错误。如何通过SCE-UA方法优化参数？

链接
顶级模型功能 | R 文档
https://www.rdocumentation.org/packages/topmodel/versions/0.7.3/topics/topmodel
R中topmodel的使用教程
https://paramo.cc.ic.ac.uk/topmodel_tutorial
topmodel：在 R
中实现水文模型 TOPMODEL https://cran.r-project.org/web/packages/topmodel/index.html
rtop 包 | R 文档
https://www.rdocumentation.org/packages/rtop/versions/0.5-14
rtop：具有可变空间支持的数据插值
https://cran.r-project.org/web/packages/rtop/index.html

install.packages("topmodel")
install.packages("Hmisc")
library(topmodel)
library(Hmisc)


# Then calculate the topographic index, the resolution should be in [m].
data(huagrahuma.dem)
DEM <- sinkfill(huagrahuma.dem, 25, 0.1)
topindex <- topidx(DEM, resolution=25)

# The values need to be split into a set of classes, since topmodel() is a semidistributed model that lumps hydrologically similar areas into the same hydrological response units.
# Here we define 16 hydrological response units:

topidx <- make.classes(topindex,16)

# the delay function is a bit more tricky because this requires cumulative fractions, but you generate it as follows:

n <- 5 # number of classes; a higher number will result in a smoother histogram
delay <- flowlength(huagrahuma.dem)*25 # TODO: add the outlet coordinates; otherwise the flowlength will be calculated to the edge of the map.
delay <- make.classes(delay, n)
delay <- delay[n:1,]
delay[,2] <- c(0, cumsum(delay[1:(n-1),2]))
############ PART 1: running the rainfall-runoff model ##############

## Load the example dataset from the Huagrahuma catchment
## and attach it to the search path

data(huagrahuma)
attach(huagrahuma)

## Initial exploration of the data:
str(huagrahuma)
topidx
parameters
rain

#I suspected the error was triggered by NA values which Qobs has. So I tried to convert NA to numerical values but ends in vain.  
k=1
for (i in Qobs) { 
  if (is.na(i)){               
    Qobs[k]<-Qobs[k-1]
  } 
  k<-k+1
}

####SCE-UA
set.seed(1)
fun=function(pars, topidx, delay, rain, ETp) topmodel:::topmodel(pars, topidx, delay, rain, ETp)
# Objective function, sud differences
OFUN = function(pars, topidx, delay, rain, ETp,Qobs) {
  Qsim=fun(pars, topidx, delay, rain, ETp)
  sum(abs(Qsim-Qobs)/Qobs)
  }
sceuares = rtop:::sceua(OFUN, pars = c(3.167914e-05,-5.990615e-01,2.129723e-02,2.626373e-02,0.05, 2.850000e+00,  1.000000e+03  ,1.199171e+03  ,9.361053e-03,7.235573e-01 ,2.500000e-01), lower = c(0,-7,0.001,0,0.005,0.001,50,50,0.0001,0,0.25), 
                        upper = c(0.0001,10,0.25,0.01,0.08,40,2000,2000,0.2,5,0.25),rain=rain,topidx=topidx,delay=delay,ETp=ETp,Qobs=Qobs)

>>>Error in if (sum(mapply(FUN = function(x, y, z) max(y - x, x - z, 0),  : 
  missing value where TRUE/FALSE needed

Answer 1

我找到了另外一个也有sceua优化功能的包，名字叫“hydromad”。通过使用它，我能够优化 TOPMODEL 参数。这是使用 huagrahuma 数据集从 topmodel 包优化 topmodel 函数时的代码。

水狂 |水文模型评估和开发
http://hydromad.catchment.org/#SCEoptim

install.packages("hydromad", repos="http://hydromad.catchment.org",dependencies = TRUE)
library(hydromad)


## Objective function
OFUN = function(pars, topidx, delay, rain, ETp,Qobs) {
  Qsim<-topmodel:::topmodel(pars, topidx, delay, rain, ETp)
  k<-1     
  cnt<-0
  obj<-0    
  for (i in Qobs) {   #In order to ignore timesteps Qobs = NA
    if (is.na(i)){               
      obj<-obj
    } else{
      obj<-obj + abs(Qsim[k] - i)/i #using relative error as Objective function
      cnt<-cnt+1
    }
    k<-k+1
  }
  obj/cnt
}

sceuares =SCEoptim(OFUN, par = c(3.167914e-05,-5.990615e-01,2.129723e-02,2.626373e-02,0.05, 2.850000e+00,  1.000000e+03  ,1.199171e+03  ,9.361053e-03,7.235573e-01 ,2.500000e-01), lower = c(0,-7,0.001,0,0.005,0.001,50,50,0.0001,0,0.25), 
         upper = c(0.0001,10,0.25,0.01,0.08,40,2000,2000,0.2,5,0.25),rain=rain,topidx=topidx,delay=delay,ETp=ETp,Qobs=Qobs, control = list())
 
 ## run the model and visualize the outcome:
 
 Qsim <- topmodel(sceuares$par, topidx, delay, rain, ETp)
 plot(Qsim, type="l", col="red")
 points(Qobs)
 
 ## Evaluate the model 
 NSeff(Qobs, Qsim) #Nash-Sutcliffe efficiency
 100*sceuares$value #MeanAbsolutePercentageError

enter image description here
带有初始参数的模型输出
enter image description here
sceua
优化参数的模型输出红线是模型输出（Qsim）。
地块是 Qobs。

Answer 2

你遇到这个问题是因为存在相同的下限和上限，正如您在脚本中定义的那样（向量的最后一个，“0.25”）：

lower = c(0,-7,0.001,0,0.005,0.001,50,50,0.0001,0,0.25)
upper = c(0.0001,10,0.25,0.01,0.08,40,2000,2000,0.2,5,0.25)

我遇到了一个类似的问题“if (sum(mapply(FUN = function(x, y, z) max(y - x, x - z, 0)”，在我改变了我的 upper 和参数的下限是正确的范围而不是相同的值，问题就解决了。

如何使用rtop包中的sceua方法优化模型参数，R语言

How to optimize model parameters using sceua method in rtop package, R language

simulation

parameters

optimization

r