hydroMOPSO: Multi-Objective Particle Swarm Optimisation algorithm (NMPSO)
In hydroMOPSO: Multi-Objective Calibration of Hydrological Models using MOPSO
View source: R/hydroMOPSOv2021.R
hydroMOPSO R Documentation
Multi-Objective Particle Swarm Optimisation algorithm (NMPSO)
Description
Multi-objective Particle Swarm Optimisation algorithm (NMPSO). The default configuration of hydroMOPSO has been adapted to obtain results with the fewest number of iterations possible.
Important: In Example 5 (calibration of GR4J hydrological model), maxit = 50 was set just for practical needs when testing the package. For acceptable results please change to maxit = 250. With any more robust model and with up to 12 parameters we recommend maxit = 1000.
Usage
hydroMOPSO(fn='hydromod',
lower=-Inf,
upper=Inf,
control=list(),
model.FUN=NULL,
model.FUN.args=list(),
obj.thr = NULL,
...)
Arguments
fn
(function or character)
Object with the name of a valid R function to be optimised (minimised or maximised). When the goal is to optimise just simple functions (problems not associated with models with input and output data), it is possible to specify the name of any function correctly defined by the user. Special cases occur when the user is working with models, declared as internal or external functions of R. In these last cases, fn='hydromod' specifies that the optimisation is applied to a model that can be invoked from R(tipically, an executable file that must be run from the system console), but is executed entirely outside of this environment. On the other hand, fn='hydromodInR' specifies that the optimisation is applied to a model that can be executed within the R environment.
In detail:
-) When fn!='hydromod' & fn!='hydromodInR', the first argument of fn has to be a vector of parameters over which optimisation is going to take place. It must return a vector with as many elements as objectives have been set in the function, and where each objective must be a scalar result. In this case, the algorithm uses the vector of values returned by fn as both model output and its corresponding set of optimised scalar results
-) When fn=='hydromod' the algorithm will optimise the R-external model defined by model.FUN and model.args, which are used to extract the values simulated by the model and to compute its corresponding goodness-of-fit measures.
-) When fn=='hydromodInR' the algorithm will optimise the R model defined by model.FUN and model.args, which are used to extract the values simulated by the model and to compute its corresponding goodness-of-fit measures.
When fn=='hydromod' | fn=='hydromodInR', the function must return a list with two (2) specific elements, the first element of the list consists of the vector with as many elements as objectives have been established in the function, and where each objective must be a scalar result; the second element of the list corresponds to a matrix with the raw output data of the model that determines the scalar results of the objectives, for example time series of a hydrological model such as streamflow, evapotranspiration, soil moisture, among others. The matrix with the raw output data of the model must have as many columns as there are simulated variables being worked on in the optimisation, and this number of variables should not necessarily coincide with the number of objectives set. for example, flows could only be returned from a hydrological model to analyze three objectives.
lower
(numeric)
Lower boundary for each parameter
In hydroMOPSO the length of lower and upper are used to defined the dimension of the solution space
upper
(numeric)
Upper boundary for each parameter
In hydroMOPSO the length of lower and upper are used to defined the dimension of the solution space
control
(list)
A list of control parameters. See‘Details’
model.FUN
(character)
(OPTIONAL) Used only when fn=='hydromod' | fn=='hydromodInR'
A valid R function representing the model code to be optimised
model.FUN.args
(list)
(OPTIONAL) Used only when fn=='hydromod' | fn=='hydromodInR'
A list with the arguments to be passed to model.FUN
obj.thr
(list)
(OPTIONAL) Used only when fn=='hydromod' | fn=='hydromodInR'
Thresholds for each objective
...
further arguments to be passed to fn
Details
By default, hydroMOPSO performs minimisation on all objectives specified in fn (MinMax='min' in control list), but this can be changed to maximisation (MinMax='max' in control list). If in fn you have to maximise some objectives and minimise others, you must make them all point to the same direction (all maximising or all minimising), which can be handled simply with a sign (See Example 2 where this type of case is presented).
Although the NMPSO algorithm was formulated to deal with many objectives, the default definitions in hydroMOPSO, and therefore the applications made in research linked to this package, have been made with a focus on two and three objectives. Extending applications to optimisations with four or more objectives is possible with this package (so you are welcome to formulate such problems and solve them with hydroMOPSO!), but be very careful in analyzing your results.
The control argument is a list that can supply any of the following components:
- drty.in
(character)
(OPTIONAL) Used only when fn='hydromod'
Name of the directory storing the input files required for PSO, i.e. ‘ParamRanges.txt’ and ‘ParamFiles.txt’.
- drty.out
(character)
Path to the directory storing the output files generated by hydroMOPSO.
- param.ranges
(character)
(OPTIONAL) Used only when fn=='hydromod' | fn=='hydromodInR'
Name of the file defining the minimum and maximum boundary values for each one of the parameters to be optimisated with NMPSO.
- digits
(numeric)
(OPTIONAL) Used only when write2disk=TRUE
Number of significant digits used for writing the output files with scientific notation.
- digits.dom
(numeric)
number of decimal places used in dominance check. Fewer decimal places (say, 16, 8, or 4, for example) may be necessary to prevent the algorithm from resulting in solutions that are nearly the same.
By default digits.dom=Inf, which basically means numbers are not rounded
- MinMax
(character)
Indicates whether a maximisation or minimisation multi-objetive problem needs to be solved. Valid values are in: c('min', 'max'). By default MinMax='min'. This control argument applies to all objective functions at the same time, so they must all go in the same direction (either all maximizing or all minimizing; keep in mind that for a particular function to go from maximizing to minimizing, or vice versa, it is only necessary add a minus sign (-)).
- npart
(numeric)
Number of particles in the swarm. By default npart=10, inherited from R-package hydroMOPSO.
- maxrep
(numeric)
Maximum number of particles to be stored in the updated Pareto Front in each iteration.
By default maxrep=100
- maxcross
(numeric)
Maximum number of Pareto Front particles that, for each iteration, perform the crossing and mutatiom in the application of genetic operators.
By default maxcross=50
- maxit
(numeric)
Maximum number of iterations.
By default maxit=1000
- Xini.type
(character)
Indicates how to initialise the particles' positions in the swarm within the ranges defined by lower and upper.
Valid values are:
-) Sobol: Sobol initialisation of positions, using npart number of samples contained in parameter space bounded by lower and upper. It requires the randtoolbox package
-) lhs: Latin Hypercube initialisation of positions, using npart number of strata to divide each parameter range. It requires the lhs package
-) random: random initialisation of positions within lower and upper
By default Xini.type='Sobol'
- Vini.type
(character)
Indicates how to initialise the particles' velocities in the swarm.
Valid values are:
-) zero: all the particles are initialised with zero velocity
-) random2011: random initialisation of velocities within lower-Xini and upper-Xini, as defined in SPSO 2011 (‘Vini=U(lower-Xini, upper-Xini)’) (see Clerc, 2012, 2010)
-) lhs2011: same as in random2011, but using a Latin Hypercube initialisation with npart number of strata instead of a random uniform distribution for each parameter. It requires the lhs package
-) random2007: random initialisation of velocities within lower and upper using the ‘half-diff’ method defined in SPSO 2007 (‘Vini=[U(lower, upper)-Xini]/2’) (see Clerc, 2012, 2010)
-) lhs2007: same as in random2007, but using a Latin Hypercube initialisation with npart number of strata instead of a random uniform distribution for each parameter. It requires the lhs package
By default Vini.type='zero'
- boundary.wall
(character)
Indicates the type of boundary condition to be applied during optimisation.
Valid values are: absorbing2011, absorbing2007, reflecting, damping, invisible
By default boundary.wall='absorbing2011'
Experience has shown that Clerc's constriction factor and the inertia weights do not always confine the particles within the solution space. To address this problem, Robinson and Rahmat-Samii (2004) and Huang and Mohan (2005) propose different boundary conditions, namely, reflecting, damping, absorbing and invisible to define how particles are treated when reaching the boundary of the searching space (see Robinson and Rahmat-Samii (2004) and Huang and Mohan (2005) for further details).
- cal.hv
(logical)
(OPTIONAL)
Indicates whether or not the hypervolume formed between the hyperplane of the Pareto Front and a nadir point designated as nadir.point will be calculated.
By default cal.hv=FALSE
- nadir.point
(numeric)
(OPTIONAL) Only required when cal.hv=TRUE
Nadir point from which the hypervolume will be calculated in each iteration step. It should correspond to a reference point considered as the worst acceptable optimal value.
- n.samples
(integer)
(OPTIONAL) Only required when cal.hv=TRUE
Number of points to estimate hypervolume, based on MonteCarlo sampling.
By default n.samples=10000
- write2disk
(logical)
Indicates if the output files will be written to the disk.
By default write2disk=TRUE
- verbose
(logical)
Indicates if progress messages are to be printed.
By default verbose=TRUE
- plot
(logical)
Indicates if a plot with the Pareto Front will be drawn after each iteration.
By default plot=FALSE
- REPORT
(integer)
(OPTIONAL) Used only when verbose=TRUE
The frequency of report messages printed to the screen.
By default REPORT=10
- parallel
(character)
Indicates how to parallelise ‘hydroMOPSO’ (to be precise, only the evaluation of the objective function fn is parallelised). Valid values are:
-)none: no parallelisation is made (this is the default value)
-)parallel: parallel computations for network clusters or machines with multiple cores or CPUs. A ‘FORK’ cluster is created with the makeForkCluster function. When fn.name='hydromod' the evaluation of the objective function fn is done with the clusterApply function of the parallel package. When fn.name != 'hydromod' the evaluation of the objective function fn is done with the parRapply function of the parallel package.
-)parallelWin: parallel computations for network clusters or machines with multiple cores or CPUs (this is the only parallel implementation that works on Windows machines). A ‘PSOCK’ cluster is created with the makeCluster function. When fn.name='hydromod' the evaluation of the objective function fn is done with the clusterApply function of the parallel package. When fn.name != 'hydromod' the evaluation of the objective function fn is done with the parRapply function of the parallel package.
- par.nnodes
(numeric)
(OPTIONAL) Used only when parallel!='none'
Indicates the number of cores/CPUs to be used in the local multi-core machine, or the number of nodes to be used in the network cluster.
By default par.nnodes is set to the amount of cores detected by the function detectCores() (parallel package)
- par.pkgs
(character)
(OPTIONAL) Used only when parallel='parallelWin'
List of package names (as characters) that need to be loaded on each node for allowing the objective function fn to be evaluated.
Value
(list)
The returned list contains elements that vary according to the input specifications
- Rep
(list)
Particle repository for the last iteration (just second last phase of NMPSO), detailing:
- Position (matrix)
Positions of each set of Pareto Front particles until the last iteration.
- Objs (matrix)
Objective values of each set of Pareto Front particles until the last iteration.
- BFE (numeric)
Balanceable Fitness Estimation (BFE) of each set of Pareto Front particles until the last iteration.
- Ranking.BFE (matrix)
Ranking of each set of Pareto Front particles until the last iteration, according to the BFE value.
- MOPSOResults
(list)
Particle repository history of all iterations (both phases of NMPSO), detailing:
- ParetoFront (data.frame)
History of objectives values of each Pareto Front particles in all iterations (both phases). In this data.frame, the first column indicates the iteration Iter; the second column the phase Phase (1 or 2); and the following columns are as many as objectives treated, being identified with the assigned name.
- Particles_ParetoFront (data.frame)
History of positions of each Pareto Front particles in all iterations (both phases). In this data.frame, the first column indicates the iteration Iter; the second column the phase Phase (1 or 2); then as many columns as objectives treated, being identified with the assigned name; and finally, as many columns as decision variables (parameters).
- MaxMin (data.frame)
Specification on whether the objectives are maximised or minimised.
- ObjsNames (data.frame)
Name of each of the objectives (Obj1, Obj2, ...).
- hydroDetails
(list)
(ONLY ADDED WHEN fn=='hydromod' | fn=='hydromodInR')
Details about the modeling involved in optimisation:
- Dimensions (data.frame)
Number of objectives and number of output variables involved in the optimisation.
- NamesAndUnitsVars (data.frame)
Name and unit of measure of the output variables involved in the optimisation (var1, var1_unit, var2, var2_unit, ...).
- Obs (list)
Observed values of each of the variables involved in the optimisation, keeping in mind that the same format indicated as mandatory input data Obs within the FUN function is maintained.
- WarmUp (data.frame)
Time series indicating the warm-up period used in the optimisation.
- DatesCal (data.frame)
Time series indicating the calibration period used in the optimisation.
- hydroResults
(list)
(ONLY ADDED WHEN fn=='hydromod' | fn=='hydromodInR')
Post-processed results about the modeling involved in optimisation:
- ParticlesFull (data.frame)
History of positions of each Pareto Front particles in all iterations. In this data.frame, the first column indicates the simulation number Sim, in ascending order from the first simulation (first iteration, phase 1) to the last simulation (last iteration, phase 2); then as many columns as objectives treated, being identified with the assigned name; and finally, as many columns as decision variables (parameters).
- FilledPOF (data.frame)
Filled Pareto front, built from evaluating the dominance of the solutions of all the iterations performed in the optimisation. To prevent the filled Pareto Front from having too many solutions, the parameters and objective values are rounded according to input DigitsDom (number of decimal places). In this data.frame, the first column indicates the simulation number Sim; then as many columns as objectives treated, being identified with the assigned name.
- ParticlesFilledPOF (data.frame)
Perticles from filled Pareto Front. In this data.frame, the first column indicates the simulation number Sim; then as many columns as objectives treated, being identified with the assigned name; and finally, as many columns as decision variables (parameters).
- ModelOut (list)
Time series of the model output variables, for all solutions of the filled Pareto Front. This list has as many objects as output variables, and each one corresponds to an object of class zoo with as many columns as solutions of the filled Pareto Front.
- ParticleBestCS (data.frame)
Best compromise solution, i.e., the solution with the minimum Euclidean distance from the maximum values of each objective. data.frame with only one row and several columns: the first column indicates the simulation number Sim; then as many columns as objectives treated, being identified with the assigned name; and finally, as many columns as decision variables (parameters).
- ModelOutBestCS (list)
Time series of the model output variables, just for the best compromise solution. This list has as many objects as output variables, and each one corresponds to an object of class zoo with a single time serie.
- ParticleBestObjs (list)
Solutions that minimise/maximise each of the objectives. data.frame with only one row. In a first level, this list has as many objects as objectives involves in the optimisation, each one with a data.frame with only one row and several columns: the first column indicates the simulation number Sim; then as many columns as objectives treated, being identified with the assigned name; and finally, as many columns as decision variables (parameters).
- ModelOutBestObjs (list)
Time series of the model output variables, for the maximisation/minimisation of each objective. In a first level, this list has as many objects as objectives involves in the optimisation and, in a second level, each one corresponds to a list with as many objects as output variables, each one corresponding to an object of class zoo with a single time serie.
- AnalysisPeriod (character)
String indicating the analysis period, in this case "calibration".
- DigitsDom (numeric)
Number of decimal places used in dominance check. Fewer decimal places (say, 16, 8, or 4, for example) may be necessary to prevent the algorithm from resulting in solutions that are nearly the same.
- ObjsNames (data.frame)
Name of each of the objectives (Obj1, Obj2, ...).
- MaxMin (data.frame)
Specification on whether the objectives are maximised or minimised, must be in c("max", "min").
- Obs (list)
Observed values of each of the variables involved in the optimisation, keeping in mind that the same format indicated as mandatory input data Obs within the FUN function is maintained.
- Dimensions (data.frame)
Number of objectives and number of output variables involved in the optimisation.
- NamesAndUnitsVars (data.frame)
Name and unit of measure of the output variables involved in the optimisation (var1, var1_unit, var2, var2_unit, ...).
- WarmUp (data.frame)
Time series indicating the warm-up period used in the optimisation.
- DatesCal (data.frame)
Time series indicating the calibration period used in the optimisation.
Note
1) For a better understanding of the application cases in which fn=='hydromod' | fn=='hydromodInR', it is strongly recommended to review the complementary tutorials to this package.
Author(s)
Rodrigo Marinao Rivas ra.marinao.rivas@gmail.com, Mauricio Zambrano-Bigiarini mzb.devel@gmail.com
References
Lin, Q., Liu, S., Zhu, Q., Tang, C., Song, R., Chen, J., Coello, C. A. C., Wong, K.-C., & Zhang, J. (2018). Particle Swarm Optimization With a Balanceable Fitness Estimation for Many-Objective Optimization Problems. IEEE Transactions on Evolutionary Computation, 22(1), 32-46. doi:10.1109/TEVC.2016.2631279
Marinao-Rivas, R., & Zambrano-Bigiarini, M. (2021). Towards best default configuration settings for NMPSO in Multiobjective Optimization. 2021 IEEE Latin American Conference on Computational Intelligence. (Accepted).
Zambrano-Bigiarini, M.; R. Rojas (2013), A model-independent Particle Swarm Optimization software for model calibration, Environmental Modelling & Software, 43, 5-25, doi:10.1016/j.envsoft.2013.01.004
Coello, C. A. C., & Lechuga, M. S. (2002). MOPSO: A proposal for multiple objective particle swarm optimization. Proceedings of the 2002 Congress on Evolutionary Computation, CEC 2002, 2, 1051-1056. doi:10.1109/CEC.2002.1004388
Kennedy, J., & Eberhart, R. (1995). Particle swarm optimization. Proceedings of ICNN'95 - International Conference on Neural Networks, 4, 1942-1948. doi:10.1109/ICNN.1995.488968
Deb, K. (1999). Multi-objective genetic algorithms: problem difficulties and construction of test problems. Evolutionary computation, 7, 205-230. doi:10.1162/EVCO.1999.7.3.205
Kursawe, F. (1991). A variant of evolution strategies for vector optimization. Lecture Notes in Computer Science (including subseries Lecture Notes in Artificial Intelligence and Lecture Notes in Bioinformatics), 496 LNCS, 193-197. doi:10.1007/BFB0029752
Deb, K., Thiele, L., Laumanns, M., & Zitzler, E. (2005). Scalable Test Problems for Evolutionary Multiobjective Optimization (bll 105-145; A. Abraham, L. Jain, & R. Goldberg, Reds). doi:10.1007/1-84628-137-7_6
Examples
###############################################################################################
# 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)
)
###############################################################################################
# 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
)
)
hydroMOPSO documentation built on April 26, 2023, 1:14 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
View source: R/hydroMOPSOv2021.R
| hydroMOPSO | R Documentation |
Multi-Objective Particle Swarm Optimisation algorithm (NMPSO)
Description
Multi-objective Particle Swarm Optimisation algorithm (NMPSO). The default configuration of hydroMOPSO has been adapted to obtain results with the fewest number of iterations possible.
Important: In Example 5 (calibration of GR4J hydrological model), maxit = 50 was set just for practical needs when testing the package. For acceptable results please change to maxit = 250. With any more robust model and with up to 12 parameters we recommend maxit = 1000.
Usage
hydroMOPSO(fn='hydromod',
lower=-Inf,
upper=Inf,
control=list(),
model.FUN=NULL,
model.FUN.args=list(),
obj.thr = NULL,
...)
Arguments
fn |
( When |
lower |
( |
upper |
( |
control |
( |
model.FUN |
( |
model.FUN.args |
( |
obj.thr |
( |
... |
further arguments to be passed to |
Details
By default, hydroMOPSO performs minimisation on all objectives specified in fn (MinMax='min' in control list), but this can be changed to maximisation (MinMax='max' in control list). If in fn you have to maximise some objectives and minimise others, you must make them all point to the same direction (all maximising or all minimising), which can be handled simply with a sign (See Example 2 where this type of case is presented).
Although the NMPSO algorithm was formulated to deal with many objectives, the default definitions in hydroMOPSO, and therefore the applications made in research linked to this package, have been made with a focus on two and three objectives. Extending applications to optimisations with four or more objectives is possible with this package (so you are welcome to formulate such problems and solve them with hydroMOPSO!), but be very careful in analyzing your results.
The control argument is a list that can supply any of the following components:
- drty.in
(
character)
(OPTIONAL) Used only whenfn='hydromod'
Name of the directory storing the input files required for PSO, i.e. ‘ParamRanges.txt’ and ‘ParamFiles.txt’.- drty.out
(
character)
Path to the directory storing the output files generated by hydroMOPSO.- param.ranges
(
character)
(OPTIONAL) Used only whenfn=='hydromod' | fn=='hydromodInR'
Name of the file defining the minimum and maximum boundary values for each one of the parameters to be optimisated with NMPSO.- digits
(
numeric)
(OPTIONAL) Used only whenwrite2disk=TRUE
Number of significant digits used for writing the output files with scientific notation.- digits.dom
(
numeric)
number of decimal places used in dominance check. Fewer decimal places (say, 16, 8, or 4, for example) may be necessary to prevent the algorithm from resulting in solutions that are nearly the same.
By defaultdigits.dom=Inf, which basically means numbers are not rounded- MinMax
(
character)
Indicates whether a maximisation or minimisation multi-objetive problem needs to be solved. Valid values are in:c('min', 'max'). By defaultMinMax='min'. Thiscontrolargument applies to all objective functions at the same time, so they must all go in the same direction (either all maximizing or all minimizing; keep in mind that for a particular function to go from maximizing to minimizing, or vice versa, it is only necessary add a minus sign (-)).- npart
(
numeric)
Number of particles in the swarm. By defaultnpart=10, inherited from R-packagehydroMOPSO.- maxrep
(
numeric)
Maximum number of particles to be stored in the updated Pareto Front in each iteration.
By defaultmaxrep=100- maxcross
(
numeric)
Maximum number of Pareto Front particles that, for each iteration, perform the crossing and mutatiom in the application of genetic operators.
By defaultmaxcross=50- maxit
(
numeric)
Maximum number of iterations.
By defaultmaxit=1000- Xini.type
(
character)
Indicates how to initialise the particles' positions in the swarm within the ranges defined bylowerandupper.
Valid values are:
-) Sobol: Sobol initialisation of positions, usingnpartnumber of samples contained in parameter space bounded bylowerandupper. It requires the randtoolbox package
-) lhs: Latin Hypercube initialisation of positions, usingnpartnumber of strata to divide each parameter range. It requires the lhs package
-) random: random initialisation of positions withinlowerandupper
By defaultXini.type='Sobol'- Vini.type
(
character)
Indicates how to initialise the particles' velocities in the swarm.
Valid values are:
-) zero: all the particles are initialised with zero velocity
-) random2011: random initialisation of velocities withinlower-Xiniandupper-Xini, as defined in SPSO 2011 (‘Vini=U(lower-Xini, upper-Xini)’) (see Clerc, 2012, 2010)
-) lhs2011: same as in random2011, but using a Latin Hypercube initialisation withnpartnumber of strata instead of a random uniform distribution for each parameter. It requires the lhs package
-) random2007: random initialisation of velocities withinlowerandupperusing the ‘half-diff’ method defined in SPSO 2007 (‘Vini=[U(lower, upper)-Xini]/2’) (see Clerc, 2012, 2010)
-) lhs2007: same as in random2007, but using a Latin Hypercube initialisation withnpartnumber of strata instead of a random uniform distribution for each parameter. It requires the lhs package
By defaultVini.type='zero'- boundary.wall
(
character)
Indicates the type of boundary condition to be applied during optimisation.
Valid values are: absorbing2011, absorbing2007, reflecting, damping, invisible
By defaultboundary.wall='absorbing2011'
Experience has shown that Clerc's constriction factor and the inertia weights do not always confine the particles within the solution space. To address this problem, Robinson and Rahmat-Samii (2004) and Huang and Mohan (2005) propose different boundary conditions, namely, reflecting, damping, absorbing and invisible to define how particles are treated when reaching the boundary of the searching space (see Robinson and Rahmat-Samii (2004) and Huang and Mohan (2005) for further details).- cal.hv
(
logical)
(OPTIONAL)
Indicates whether or not the hypervolume formed between the hyperplane of the Pareto Front and a nadir point designated asnadir.pointwill be calculated.
By defaultcal.hv=FALSE- nadir.point
(
numeric)
(OPTIONAL) Only required whencal.hv=TRUE
Nadir point from which the hypervolume will be calculated in each iteration step. It should correspond to a reference point considered as the worst acceptable optimal value.- n.samples
(
integer)
(OPTIONAL) Only required whencal.hv=TRUE
Number of points to estimate hypervolume, based on MonteCarlo sampling.
By defaultn.samples=10000- write2disk
(
logical)
Indicates if the output files will be written to the disk.
By defaultwrite2disk=TRUE- verbose
(
logical)
Indicates if progress messages are to be printed.
By defaultverbose=TRUE- plot
(
logical)
Indicates if a plot with the Pareto Front will be drawn after each iteration.
By defaultplot=FALSE- REPORT
(
integer)
(OPTIONAL) Used only whenverbose=TRUE
The frequency of report messages printed to the screen.
By defaultREPORT=10- parallel
(
character)
Indicates how to parallelise ‘hydroMOPSO’ (to be precise, only the evaluation of the objective functionfnis parallelised). Valid values are:
-)none: no parallelisation is made (this is the default value)
-)parallel: parallel computations for network clusters or machines with multiple cores or CPUs. A ‘FORK’ cluster is created with themakeForkClusterfunction. Whenfn.name='hydromod'the evaluation of the objective functionfnis done with theclusterApplyfunction of the parallel package. Whenfn.name != 'hydromod'the evaluation of the objective functionfnis done with theparRapplyfunction of the parallel package.
-)parallelWin: parallel computations for network clusters or machines with multiple cores or CPUs (this is the only parallel implementation that works on Windows machines). A ‘PSOCK’ cluster is created with themakeClusterfunction. Whenfn.name='hydromod'the evaluation of the objective functionfnis done with theclusterApplyfunction of the parallel package. Whenfn.name != 'hydromod'the evaluation of the objective functionfnis done with theparRapplyfunction of the parallel package.- par.nnodes
(
numeric)
(OPTIONAL) Used only whenparallel!='none'
Indicates the number of cores/CPUs to be used in the local multi-core machine, or the number of nodes to be used in the network cluster.
By defaultpar.nnodesis set to the amount of cores detected by the functiondetectCores()(parallel package)- par.pkgs
(
character)
(OPTIONAL) Used only whenparallel='parallelWin'
List of package names (as characters) that need to be loaded on each node for allowing the objective functionfnto be evaluated.
Value
(list)
The returned list contains elements that vary according to the input specifications
- Rep
(
list)
Particle repository for the last iteration (just second last phase of NMPSO), detailing:
- Position (
matrix)
Positions of each set of Pareto Front particles until the last iteration.
- Objs (matrix)
Objective values of each set of Pareto Front particles until the last iteration.
- BFE (numeric)
Balanceable Fitness Estimation (BFE) of each set of Pareto Front particles until the last iteration.
- Ranking.BFE (matrix)
Ranking of each set of Pareto Front particles until the last iteration, according to the BFE value.
- MOPSOResults
(
list)
Particle repository history of all iterations (both phases of NMPSO), detailing:
- ParetoFront (
data.frame)
History of objectives values of each Pareto Front particles in all iterations (both phases). In thisdata.frame, the first column indicates the iterationIter; the second column the phasePhase(1 or 2); and the following columns are as many as objectives treated, being identified with the assigned name.
- Particles_ParetoFront (data.frame)
History of positions of each Pareto Front particles in all iterations (both phases). In thisdata.frame, the first column indicates the iterationIter; the second column the phasePhase(1 or 2); then as many columns as objectives treated, being identified with the assigned name; and finally, as many columns as decision variables (parameters).
- MaxMin (data.frame)
Specification on whether the objectives are maximised or minimised.
- ObjsNames (data.frame)
Name of each of the objectives (Obj1, Obj2, ...).
- hydroDetails
(
list)
(ONLY ADDED WHENfn=='hydromod' | fn=='hydromodInR')
Details about the modeling involved in optimisation:
- Dimensions (
data.frame)
Number of objectives and number of output variables involved in the optimisation.
- NamesAndUnitsVars (data.frame)
Name and unit of measure of the output variables involved in the optimisation (var1, var1_unit, var2, var2_unit, ...).
- Obs (list)
Observed values of each of the variables involved in the optimisation, keeping in mind that the same format indicated as mandatory input dataObswithin theFUNfunction is maintained.
- WarmUp (data.frame)
Time series indicating the warm-up period used in the optimisation.
- DatesCal (data.frame)
Time series indicating the calibration period used in the optimisation.
- hydroResults
(
list)
(ONLY ADDED WHENfn=='hydromod' | fn=='hydromodInR')
Post-processed results about the modeling involved in optimisation:
- ParticlesFull (
data.frame)
History of positions of each Pareto Front particles in all iterations. In thisdata.frame, the first column indicates the simulation numberSim, in ascending order from the first simulation (first iteration, phase 1) to the last simulation (last iteration, phase 2); then as many columns as objectives treated, being identified with the assigned name; and finally, as many columns as decision variables (parameters).
- FilledPOF (data.frame)
Filled Pareto front, built from evaluating the dominance of the solutions of all the iterations performed in the optimisation. To prevent the filled Pareto Front from having too many solutions, the parameters and objective values are rounded according to inputDigitsDom(number of decimal places). In thisdata.frame, the first column indicates the simulation numberSim; then as many columns as objectives treated, being identified with the assigned name.
- ParticlesFilledPOF (data.frame)
Perticles from filled Pareto Front. In thisdata.frame, the first column indicates the simulation numberSim; then as many columns as objectives treated, being identified with the assigned name; and finally, as many columns as decision variables (parameters).
- ModelOut (list)
Time series of the model output variables, for all solutions of the filled Pareto Front. This list has as many objects as output variables, and each one corresponds to an object of class zoo with as many columns as solutions of the filled Pareto Front.
- ParticleBestCS (data.frame)
Best compromise solution, i.e., the solution with the minimum Euclidean distance from the maximum values of each objective. data.frame with only one row and several columns: the first column indicates the simulation numberSim; then as many columns as objectives treated, being identified with the assigned name; and finally, as many columns as decision variables (parameters).
- ModelOutBestCS (list)
Time series of the model output variables, just for the best compromise solution. This list has as many objects as output variables, and each one corresponds to an object of class zoo with a single time serie.
- ParticleBestObjs (list)
Solutions that minimise/maximise each of the objectives. data.frame with only one row. In a first level, this list has as many objects as objectives involves in the optimisation, each one with a data.frame with only one row and several columns: the first column indicates the simulation numberSim; then as many columns as objectives treated, being identified with the assigned name; and finally, as many columns as decision variables (parameters).
- ModelOutBestObjs (list)
Time series of the model output variables, for the maximisation/minimisation of each objective. In a first level, this list has as many objects as objectives involves in the optimisation and, in a second level, each one corresponds to a list with as many objects as output variables, each one corresponding to an object of class zoo with a single time serie.
- AnalysisPeriod (character)
String indicating the analysis period, in this case"calibration".
- DigitsDom (numeric)
Number of decimal places used in dominance check. Fewer decimal places (say, 16, 8, or 4, for example) may be necessary to prevent the algorithm from resulting in solutions that are nearly the same.
- ObjsNames (data.frame)
Name of each of the objectives (Obj1, Obj2, ...).
- MaxMin (data.frame)
Specification on whether the objectives are maximised or minimised, must be inc("max", "min").
- Obs (list)
Observed values of each of the variables involved in the optimisation, keeping in mind that the same format indicated as mandatory input dataObswithin theFUNfunction is maintained.
- Dimensions (data.frame)
Number of objectives and number of output variables involved in the optimisation.
- NamesAndUnitsVars (data.frame)
Name and unit of measure of the output variables involved in the optimisation (var1, var1_unit, var2, var2_unit, ...).
- WarmUp (data.frame)
Time series indicating the warm-up period used in the optimisation.
- DatesCal (data.frame)
Time series indicating the calibration period used in the optimisation.
Note
1) For a better understanding of the application cases in which fn=='hydromod' | fn=='hydromodInR', it is strongly recommended to review the complementary tutorials to this package.
Author(s)
Rodrigo Marinao Rivas ra.marinao.rivas@gmail.com, Mauricio Zambrano-Bigiarini mzb.devel@gmail.com
References
Lin, Q., Liu, S., Zhu, Q., Tang, C., Song, R., Chen, J., Coello, C. A. C., Wong, K.-C., & Zhang, J. (2018). Particle Swarm Optimization With a Balanceable Fitness Estimation for Many-Objective Optimization Problems. IEEE Transactions on Evolutionary Computation, 22(1), 32-46. doi:10.1109/TEVC.2016.2631279
Marinao-Rivas, R., & Zambrano-Bigiarini, M. (2021). Towards best default configuration settings for NMPSO in Multiobjective Optimization. 2021 IEEE Latin American Conference on Computational Intelligence. (Accepted).
Zambrano-Bigiarini, M.; R. Rojas (2013), A model-independent Particle Swarm Optimization software for model calibration, Environmental Modelling & Software, 43, 5-25, doi:10.1016/j.envsoft.2013.01.004
Coello, C. A. C., & Lechuga, M. S. (2002). MOPSO: A proposal for multiple objective particle swarm optimization. Proceedings of the 2002 Congress on Evolutionary Computation, CEC 2002, 2, 1051-1056. doi:10.1109/CEC.2002.1004388
Kennedy, J., & Eberhart, R. (1995). Particle swarm optimization. Proceedings of ICNN'95 - International Conference on Neural Networks, 4, 1942-1948. doi:10.1109/ICNN.1995.488968
Deb, K. (1999). Multi-objective genetic algorithms: problem difficulties and construction of test problems. Evolutionary computation, 7, 205-230. doi:10.1162/EVCO.1999.7.3.205
Kursawe, F. (1991). A variant of evolution strategies for vector optimization. Lecture Notes in Computer Science (including subseries Lecture Notes in Artificial Intelligence and Lecture Notes in Bioinformatics), 496 LNCS, 193-197. doi:10.1007/BFB0029752
Deb, K., Thiele, L., Laumanns, M., & Zitzler, E. (2005). Scalable Test Problems for Evolutionary Multiobjective Optimization (bll 105-145; A. Abraham, L. Jain, & R. Goldberg, Reds). doi:10.1007/1-84628-137-7_6
Examples
###############################################################################################
# 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)
)
###############################################################################################
# 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
)
)
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.