###############################################################################################
# Example 1. Basic Benchmark function in minimisation
###############################################################################################
# This basic Benchmark function has two objectives (M = 2) in minimisation, its Pareto optimal
# front is discontinuous with four disconnected curves.This function works with 2 decision
# variables (D = 2).
# Main reference for function: Deb (1999)
library(hydroMOPSO)
lower <- c(0, 0)
upper <- c(1, 1)
fnBasic <- function(param){
x1 <- param[1]
x2 <- param[2]
obj1 <- x1
obj2 <- (1 + 10*x2)*(1-(x1/(1+10*x2))^2 - x1/(1+10*x2)*sin(2*pi*4*x1))
out <- list(c(obj1, obj2)) # For consistency with further examples, this must be a list
names(out) <- "Objs" # The name "Objs" is a mandatory requirement
return(out)
}
set.seed(100) # Setting the seed (for reproducible results)
out <- hydroMOPSO(fn = fnBasic,
lower = lower,
upper = upper,
control=list(npart = 10, maxrep = 100, maxcross = 50,
MinMax = "min", maxit = 50, plot = TRUE)
)
# \donttest{
###############################################################################################
# Example 2. Basic Benchmark function in maximisation
###############################################################################################
#
# This example is identical to Example 1, but the functions are in maximisation
#
# IMPORTANT:
# In the literature related to multi-objective optimisation, test functions are usually
# presented with minimisation objectives (such as the function used in Example 1). However,
# this does not necessarily always have to be formulated that way, especially when it comes to
# real-world applications.
#
# In this second example we just want to remind you that the disjunctive between maximising or
# minimising objectives is "a matter of signs".
#
# With this in account, as explained in the documentation, for hydroMOPSO operation there is
# one requirement which you MUST TAKE CARE OF:
#
# "The problems must be formulated in such a way that ALL objectives are either maximising or
# minimising. If your problem mixes both types of objectives, just add minus signs (-) in the
# results that require it..."
#
# Main reference for function: Kursawe (1991)
library(hydroMOPSO)
lower <- c(0, 0)
upper <- c(1, 1)
fnBasic <- function(param){
x1 <- param[1]
x2 <- param[2]
obj1 <- -( x1 )
obj2 <- -( (1 + 10*x2)*(1-(x1/(1+10*x2))^2 - x1/(1+10*x2)*sin(2*pi*4*x1)) )
# note tha minus sign was added in obj1 and obj2
out <- list(c(obj1, obj2)) # For consistency with further examples, this must be a list
names(out) <- "Objs" # The name "Objs" is a mandatory requirement
return(out)
}
set.seed(100) # Setting the seed (for reproducible results)
out <- hydroMOPSO(fn = fnBasic,
lower = lower,
upper = upper,
control=list(npart = 10, maxrep = 100, maxcross = 50,
MinMax = "max", maxit = 50, plot = TRUE) # note that now MinMax="max"
)
###############################################################################################
# Example 3. Using 'smoof' package: Kursawe function
###############################################################################################
#
# This Benchmark function has two objectives (M = 2), its Pareto optimal front is discontinuous
# and non-convex. For this example it will be implemented with 3 decision variables (D = 3)
#
# Main reference for function: Kursawe (1991)
library(hydroMOPSO)
library(smoof)
D <- 3
lower <- rep(-5,D)
upper <- rep(5,D)
Kursawe <- smoof::makeKursaweFunction(D) # using 'smoof' package
fnKursawe <- function(param){
objs <- Kursawe(x = param)
obj1 <- objs[1]
obj2 <- objs[2]
out <- list(c(obj1, obj2)) # For consistency with further examples, this must be a list
names(out) <- "Objs" # The name "Objs" is a mandatory requirement
return(out)
}
set.seed(100) # Setting the seed (for reproducible results)
out <- hydroMOPSO(fn = fnKursawe,
lower = lower,
upper = upper,
control=list(npart = 10, maxrep = 100, maxcross = 50,
MinMax = "min", maxit = 50, plot = TRUE)
)
###############################################################################################
# Example 4. Using 'smoof' package: DTLZ2 function with three objectives
###############################################################################################
#
# In this example, this Benchmark is formulated with two objectives (M = 3) and 12 decision
# variables (D = 12) its Pareto optimal front is concave.
#
# Main reference for function: Deb (2005)
library(hydroMOPSO)
library(smoof)
M <- 3
D <- 12
lower <- rep(0,D)
upper <- rep(1,D)
DTLZ2 <- smoof::makeDTLZ2Function(D, M) # using 'smoof' package
fnDTLZ2 <- function(param){
objs <- DTLZ2(x = param)
obj1 <- objs[1]
obj2 <- objs[2]
obj3 <- objs[3]
out <- list(c(obj1, obj2, obj3)) # For consistency with further examples, this must be a list
names(out) <- "Objs" # The name "Objs" is a mandatory requirement
return(out)
}
set.seed(100) # Setting the seed (for reproducible results)
out <- hydroMOPSO(fn = fnDTLZ2,
lower = lower,
upper = upper,
control=list(npart = 10, maxrep = 100, maxcross = 50,
MinMax = "min", maxit = 50, plot = TRUE)
)
###############################################################################################
# Example 5. Calibration of GR4J hydrological model
###############################################################################################
#
# For this example, a "real-world" problem has been formulated: the calibration of a
# hydrological model
#
# In detail...
# Hydrological model: GR4J (Perrin et al., 2004)
# Number of parameters: four (X1, X2, X3, X4; see Perrin et al. (2004))
# Study area: Trancura River Basin (RTL)
# Input variables: Precipitation (pcp) and Potetntial EvapoTranspiration (pet)
# Calibration output variable: Streamflow (qobs)
library(hydroMOPSO)
library(airGR)
library(hydroTSM)
library(hydroGOF)
# RTL basin ------------------------------------------------
basin.area <- 1415025887 # basin area in square meters
# Load time series -----------------------------------------
data(Trancura9414001plus) # Load RTL data set
# Dates ----------------------------------------------------
dates.raw <- Trancura9414001plus[,"Date"]
dates <- as.Date(dates.raw) # dates
# INPUTS time series ---------------------------------------
# Precipitation (input variable)
ts.pcp.raw <- Trancura9414001plus[,"P_mm"]
ts.pcp <- zoo(ts.pcp.raw, dates)
# Potential EvapoTranspiration (input variable)
ts.pet.raw <- Trancura9414001plus[,"PET_mm"]
ts.pet <- zoo(ts.pet.raw, dates)
# OUTPUTS time series --------------------------------------
# Observed streamflow (calibration output variable)
ts.qobs.raw <- Trancura9414001plus[,"Qobs_m3s"]
ts.qobs <- zoo(ts.qobs.raw, dates)
# Parameter ranges and noCal parameters --------------------
lower <- c("X1" = 0.01, "X2" = -100, "X3" = 0.01, "X4" = 0.5) # parameter range lower threshold
upper <- c("X1" = 1200, "X3" = 100, "X3" = 5000, "X4" = 5 ) # parameter range upper threshold
noCal.param <- (lower + upper)/2 # uncalibrated parameters
# Names and units of observed output variables -------------
char.obs.names <- "Streamflow"
char.obs.units <- "m3/s"
# Objectives names -----------------------------------------
char.objs.names <- c("KGE2012_Q", "KGEGarcia_Q")
# Calibration dates and subsetting -------------------------
WarmUpCal.dates <- dip("1979-01-01", "1979-12-31") # WarmUp for Calibration
Cal.dates <- dip("1980-01-01", "1999-12-31") # Calibration
FullCal.dates <- dip("1979-01-01", "1999-12-31") # WarmUp + Calibration
start.FullCal <- FullCal.dates[1]
end.FullCal <- FullCal.dates[length(FullCal.dates)]
# INPUTS time series ---------------------------------------
# Precipitation (input variable)
ts.pcp.FullCal <- window(ts.pcp, start = start.FullCal, end = end.FullCal) # subsetting pcp
# Potential EvapoTranspiration (input variable)
ts.pet.FullCal <- window(ts.pet, start = start.FullCal, end = end.FullCal) # subsetting pet
# OUTPUTS time series --------------------------------------
# Observed streamflow (calibration output variable)
ts.qobs.FullCal <- window(ts.qobs, start = start.FullCal, end = end.FullCal) # subsetting qobs
list.obs.Cal <- list(Q = ts.qobs.FullCal)
# Structuring Inputs and Options of GR4J model -------------
InputsModel.Cal <- CreateInputsModel(FUN_MOD= RunModel_GR4J,
DatesR= as.POSIXlt(FullCal.dates),
Precip= coredata(ts.pcp.FullCal),
PotEvap= coredata(ts.pet.FullCal))
RunOptions.Cal <- CreateRunOptions(FUN_MOD= RunModel_GR4J, InputsModel= InputsModel.Cal,
IndPeriod_Run = 1:length(FullCal.dates), warnings = FALSE)
# hydroMOPSO calibration -----------------------------------
set.seed(100) # Setting the seed (for reproducible results)
Cal.results <- hydroMOPSO(fn="hydromodInR",
lower=lower,
upper=upper,
control=list(MinMax="max", Xini.type = "lhs", npart=10,
maxit=50, # for better results set maxit=250 in this case
maxrep = 100, maxcross = 50,
maxeval = 15000, write2disk = FALSE, REPORT=1,
digits = 8, plot = TRUE, parallel = "none"),
model.FUN="GR4JExampleCal",
model.FUN.args = list(Obs=list.obs.Cal, # mandatory
Objs.names = char.objs.names, # mandatory
var.names = char.obs.names, # mandatory
var.units = char.obs.units, # mandatory
full.period = FullCal.dates, # mandatory
warmup.period = WarmUpCal.dates,
cal.period = Cal.dates,
# Model specific inputs
InputsModel = InputsModel.Cal, # model specific
RunOptions = RunOptions.Cal, # model specific
area = basin.area # model specific
)
)
# }
Run the code above in your browser using DataLab