[R-sig-dyn-mod] Qn about FME package
Thomas Petzoldt
thomas.petzoldt at tu-dresden.de
Wed Jul 20 21:12:42 CEST 2016
I guess the error may be in:
Sfun=sensFun(gwcost,parms=pars,sensvar=c("h"))
See example of ?sensFun, where func is a wrapper around the model that
returns the modelled (i.e. dependent or state) variables. In the
statement above, it was called with a cost function.
Thomas P.
Am 20.07.2016 um 18:02 schrieb Vishal Mehta:
> #########################################
> #Jan 2016: Dr. Satkumar Tomer's 1-D model copied
> ####GW model 1-D GW model function
> gw_model <- function(pars, hmin.under.flow, pd.under.flow, rain){
> n <- length(rain)
> net_recharge <- rf*rain
> h <- array(0,n)
> discharge <- array(0,n)
> time = array(0,n)
> h[1] <- hini
> time[1] = 1
> for (i in 1:(n-1)){
> discharge[i] <- sy*(1-pd.under.flow)*(h[i]-hmin.under.flow+ net_recharge[i]/sy)
> if (discharge[i]<0){
> discharge[i] <- 0
> }
> h[i+1] <- h[i] + (net_recharge[i]-discharge[i])/sy
> time[i+1] = time[i] + 1
> }
> return(data.frame(time=time,h=h,discharge=discharge))
> }
> #synthetic forcing
> n = 12*20 # temporal length in months
> set.seed(123) # to have the repeatibility in results
> rain <- rgamma(n, shape=8, rate = 0.1) # synthetic rainfall
> #barplot(rain,xlab="Time",ylab="Rainfall (mm)",main="Monthly rainfall")
> rain=rain/1000#mm to m of rain
> #new
>
>
> #
> #input parameters
> rf <- 0.1 # recharge factor, fraction
> sy <- 0.005 #specific yield, dimensionless
> hini <- 150 # initial condition, metres
> hmin.under.flow <- 100
> pd.under.flow<-0.967
> pars=c("hini"=hini,"rf"=rf,"sy"=sy)##this is the subset of parameters to be SA'd, calibrated, and to inform UA
> ##
> pdf("output.pdf")
> output=gw_model(pars,hmin.under.flow,pd.under.flow,rain)
> plot(output$h,type='l',xlab='time',ylab='simulated gw level(m)')
> dev.off()
> ##########################
> #July 15, 2016
> #trial following FME tutorial
> #########
> library(FME)
> #Create dummy observations
> #n=240 length of simulated output
> #output$h contains the simulated head
> #below, observation dataframe is called DataT. It has to(?) contain a "time" and "h" variable like the simulation output
> DataT=data.frame(cbind("time"=seq(1,length(n)),"h"=output$h+rnorm(sd=0.5,mean=5,n=240),"sd"=4.5))
> DataT2=DataT[,c(1,2)]
> ##
> #Create a function to calculate fit
> gwcost=function(pars){
> out=gw_model(pars,hmin.under.flow, pd.under.flow, rain)
> #cost=modCost(model=out,obs=DataT,err="sd")
> cost=modCost(model=out,obs=DataT2)
> return(cost=cost)
> }
> ###Above seems to be working(maybe not quite?)
> gwcost(pars)
> gwcost(pars)$model
> #[1] 1075 with no weighting
> plot(gwcost(pars)$residuals$res)
> #
> ##Sensitivity using sensFun (here onwards results are not following the JSS article)
> ##local sensitivity does not seem to be working. does it work only with differential equations?
> Sfun=sensFun(gwcost,parms=pars,sensvar=c("h"))
>
> #parameter tinyy to add perturbation? not making any difference
> summary(Sfun)
> plot(Sfun)
> #########
More information about the R-sig-dynamic-models
mailing list