My very first try to optimize APSIM parameters with R
Hello all,
It's me again with a short tutorial on how to use apsimr package to interact with APSIM model and use optim function in R to optimize some of the APSIM input variables. For those of you, who are not familiar with APSIM model, here is the official definition:
"The Agricultural Production Systems sIMulator (APSIM) is internationally recognised as a highly advanced simulator of agricultural systems. It contains a suite of modules which enable the simulation of systems that cover a range of plant, animal, soil, climate and management interactions." (https://www.apsim.info).
Actually, I attended in set of APSIM related presentations last week in our department, which a couple of users seemed to be trying to conceptually make their model work better (whatever we wanna call it, either spining up, calibration or optimization). Even though they seemed to know how to use 'apsimr' package to interact with the APSIM model through R, and they knew how to run some sensitivity or uncertainty analysis, I didn't see any of them trying to perform some optimizations.
For those of you also who doesn't know about 'apsimr', it's package developed by bryan stanfill (https://github.com/stanfill) with the aim of running apsim simulations through R.
So, long story short in this package you have access to the two fundamental functions for running APSIM simulations. These functions are 'apsim' and 'edit_apsim' functions that they run and edit apsim files respectively. But bryan has also provided a vectorized version of these two functions which it lets you to do both simultaneously with a just one function and it's called 'apsim_vector'. This actually lets the user to plug in this function as a model to the sensitivity analysis functions, like soboljansen function from sensitivity package. Like this :
#############################################################
"soboljansen(model = apsim_vector, X1, X2, exe = apsimExe, wd = apsimWd, vars = vars, to.run = apsimFile, g = meanYield, overwrite = TRUE)"
##################################################################
So I used this capability to run 'optim' function. In the follwoing code, I tried to optimize SummerConna and SummerU parameters by minimizing an index which is sum of squers of the difference between observed (made up values) and simulated yield. Here is the code:
I just spent less than half an hour to write this code and basically it's my first try in this matter. In the objective function, I simply tried to feed the the 'apsim_vector' wih the new parameters sent by 'optim' function. Then I compared the model's output with some made up obeservation by making an index like sums of squers which needs to be minimized.
I guess that's it !
Please share with me any comments or any idea that could help to improve the code
and happy charismasin advance
######################################################## Libraries library(apsimr) #################### Setting up the eniviroment for apsim r apsimExe <-"C:/Program Files (x86)/Apsim77-r3615/Model/Apsim.exe" apsimWd <- "~/APSIM" apsimmod<-"C:/Program Files (x86)/Apsim77-r3615/Model/" apsimFile <- "Centro.apsim" obs<-c(900,1200) #########################################objective function # ## This function just take out the yield from the output meanCowpea<-function(X){ return((X$yield)) }
######## Objective function objfun <- function(prms, odata,apsimVar){ ## I need this to run the apsim for every try apsimValue<-matrix(sapply(prms,function(x){rep(x,length(apsimVar))}),nrow=length(apsimVar),byrow=T) uniRes <- apsim_vector(X = apsimValue, exe = apsimExe, wd = apsimWd, vars = apsimVar, to.run = apsimFile, to.edit = apsimFile, g = meanCowpea) ### Arbitary index rss<-sum(((uniRes)-odata)^2) rss }
############################################# ap.var<- c("SoilWater/SummerCona", "SoilWater/SummerU") ### just to check and see how the function works objfun(c(2,3),odata = obs,ap.var)
### real optim function op <- optim(c(2,3), objfun, odata = obs,apsimVar=ap.var,control = list(trace = 1))
########################################################################
