VIC model Optimization

This page describes the Multi-Objective optimization method for calibration the VIC model.  Before conducting an automated calibration, users should familiarize themselves with standard VIC model calibration procedures.  Automated calibration is a tool that can help with obtaining better simulation results, but is not a replacement for understanding how the model responds to its parameters, and understanding the physical and theoretical meaning of those parameters. 

Multi-Objective Optimization

Multi-objective optimization moves beyond standard optimization methods that try to minimize a single objective (say Nash-Sutcliffe Efficiency) by providing a methodology for optimizing several objective functions.  Single objective optimization is limited by the choice of a single objective  function.  Every objective function will have its own set of sensitivities and limitations.  Returning for a moment to the Nash-Sutcliffe Efficiency, this objective function tends to be more sensitive to how well the model replicates observed extremes (especially peak flows), then about the models ability to reproduce the volume of observed streamflow.  Therefore a good value of Nash-Sutcliffe Efficiency could be achieved by calibrating to match peak flows, while simulating half or perhaps twice as much flow as is observed on an annual basis.  A multi-objective optimization method would let the user optimize the model for both the Nash-Sutcliffe Efficiency and Mean Annual Bias so that a set of parameters is identified that balances capturing the peak with getting mean annual flow volumes about right.  Multi-objective optimization methods can also make use of objective functions for other variables, say to optimize snow depth and streamflow, or streamflow and nutrient transport.  As long as there are observational data sets for variables simulated by the model, and an objective function can be defined, the multi-objective optimization method can be applied.

The biggest drawback of the multi-objective method of optimization is that it will not result in a single set of parameters that will minimize the objective function.  Instead they will generally result in a population of optimized parameters that will form what is called a Pareto set.  In the Pareto set, each set of parameters provides a slightly different trade-off between minimizing the various objective functions, but no set of parameters optimizes the set of objective functions better than any other. An example showing the difference between parameter space and objective space as taken from Yapo et al. (1998) is shown below.  Here the parameter X1 is allowed to vary from 0 (point B) to 1 (point A), while the parameter X2 is locked at a value of 0.  In objective space, as the value of X1 changes from A to B it shifts from minimizing objective function F2 to minimizing objective function F1.  Values of X1 between A and B populate the Pareto where neither objective function is optimized independently, but they are optimized jointly.

Here is a link to the source paper for the MOCOM-UA algorithm

Requirements for MOCOM-UA Optimization

The optimization repository contains the programs and methods required to apply the multi-objective optimization routine from the University of Arizona (MOCOM-UA) to VIC model watersheds.  The MOCOM-UA program is a general wrapper, so control scripts can be written to optimize almost anything assuming that it can handle the list of new parameters from the MOCOM-UA program and return the values of the multi-objective functions.  The control scripts provided as part of the GitHub repository are designed to facilitate running full VIC model setups, including changing VIC model soil and lake parameter files to reflect changes to input parameters, running the VIC model, running the routing model and computing objective functions relating simulated to observed streamflow. 

Accessing the Optimizer and Supporting Code

The MOCOM-UA source code and associated files are available from the PHIG GitHUB site, https://github.itap.purdue.edu/PHIGOrganization/Optimization.  The latest MOCOM-UA code is compiled and available in /depot/phig/apps/linux for PHIG team members with access to the Data Depot and Negishi cluster accounts.  

Getting Started with the Optimizer

The following instructions assume that you are already comfortable with and have a working setup of the VIC model, and that you are now ready to make use of the optimizer to improve your simulation parameterization.  While code development has been focused on minimizing the effort needed to get the optimizer set up, this is by no means an automated process.

In general the process is to:

• Establish a base simulation, so a folder that only contains what is needed to run and analyze the VIC model and to control the optimization process.  Files that must change during calibration/optimization go here (typically the global file, soil parameter file and the lake parameter file if being used).  Other files, such as the vegetation parameter and library files, can be stored in this folder but will be copied for each calibration/optimization simulation.  Large files should be stored on scratch (for speed) and referred to as absolute paths, otherwise copying these files slows down the entire optimization process unnecessarily.

• Write a simulation control script that will (1) run the VIC model using parameters and files in the base simulation; (2) post process simulation output, for example apply the routing model; and (3) calculate and report the results of objective functions that compare the current simulation to select files of "observations".

• Run the calibrator or optimizer to run the model for multiple parameter sets.  The calibrator prompts the user for each new set of parameters and runs the model using the simulation control script for the selected parameters.  The optimizer automates the selection of parameters and uses the MOCOM process to find the optimal set of parameter sets (the Pareto set).  The base simulation should be copied to scratch before starting the calibrator or optimizer.  Both tools will build a complex structure of sub directories representing each new generation or calibration parameter set and each member of that generation (for the optimizer).  The optimizer will keep track of all parameter sets and objective function results in the optimization log file, but only keep a subset of the simulation results - you control how many of the intermediate generations are saved.  Each saved folder will contain a submission script that can be used to rerun the VIC model simulation for that folder. 

Here are much more detailed steps for getting started:

Step 1: Establishing the base simulation

Within your Home Drive or Data Depot, create a new directory that will serve as the base simulation for the optimization process. You want this in a secure, backed up location so that all setup customizations are preserved., but will copy the contents to the Scratch Drive before running the calibration or optimization process.

Basic VIC model setup

The base simulation folder should include all files needed to successfully run the VIC model for you desired calibration process. When deciding which files to copy to this directory consider the following:

• Any file that must be changed in the calibration process MUST be in this folder, typically this includes the global parameter file, the soil parameter file and lake and wetland parameter file if being used.

• It is better to put large files or large numbers of smaller files that will not change in a fixed location on the relevant Scratch Drive. These can then be referred to using absolute paths, and will not have to be copied for each new set of parameters, which can significantly slow down your process.

• Other files can be stored in the base simulation folder, or in a fixed location - it is your choice.

• Be careful about using absolute paths to files stored on the Depot Drive or your Home Drive as these connections are slower than the Scratch Drive. Additionally, recognize that files on the Scratch Drive are purged regularly, so every file located on the Scratch Drive will need to be reconstituted before the optimizer can be reused after a hiatus.

Edit the VIC model global control file to reflect any new locations and start the model from the command line to confirm that it starts. Results from the VIC model simulation should be stored in the current directory (e.g., ./). If you want to test the entire simulation, open an interactive shell on the cluster and run the code there - not on the front end machine.

Post processing setup

Next, if you are going to run the routing model or another application on the VIC model results as part of the assessment process, you need to make sure that everything required to complete that analysis is available from the base scenario folder. Again, consider the size of number of files when deciding whether to copy relevant files to the base scenario folder. Local files can be used with relative paths, while those in fixed locations should use absolute paths.

Make sure that your post processing tools work from the base scenario folder, and that output from these steps is stored in the current directory.

Optimizer setup

Finally, you will want to copy the contents of the run_scripts folder from the optimizer GitHub repository into the base scenario folder. Best way to do this is to clone the repository to your personal GitHub folder, then copy the contents of the run_scripts directory to base scenario folder. Then you can keep the cloned repository up to date with GitHub, and maintain your customized code in a separate directory.

Step 2: Create a script to compute simulation statistics (objective function metrics)

With Step 1 complete, you should have all of your simulation and post processing output files stored locally in the base scenario folder. The next step is to write a script that will use the output from the model and post processor and compute values for each of the objective functions selected to evaluate the current parameter set. As this is a multi-objective optimization process, typically you will want 2 or more objective functions (every additional objective function will increase the number of simulations required to find a solution, so not too many). For example, for optimization versus river discharge, I often use Nash-Sutcliffe Efficiency (NSE), Peak Difference and the Number of Sign Counts or mean absolute bias. Three objective functions the measure differences between the simulated and observed datasets using different characteristics. Avoid using multiple objective functions that measure the same characteristics, such as Nash-Sutcliffe Efficiency and Correlation Coefficient (R-square).

The code CalculateSimulationStat.py provided with the Optimization run scripts is a good starting point, especially if you are going to optimize discharge since it is already configured to strip out matching periods and apply three objective functions. It also reads critical information about the current situation environment from the optimization control file (Step 3), so that it can be used with a variety of model setups with minimal change. If you do write your own function, consider doing the same for increased flexibility.

No matter what script you use, it must do the following:

• Process VIC model and post process output from the current directory.

• Compute objective functions that are minimized as the simulated and observed datasets get closer to each other.

  • This may mean that you need to modify the objective function, for example NSE falls between 1 and negative infinity, with 1 indicating perfect agreement between the simulated and observed datasets. The optimizer tries to minimize the objective function, which in this case would mean making the model perform as poorly as possible. CalculateSimulationStat.py works around this by reporting the value 1 - NSE, so 0 is a perfect simulation and positive infinity is a terrible simulation.

  • Also watch out for objective functions where the best value is zero, such a bias. The optimizer is minimizing the value, not finding the root (or zero value), so it will preferentially seek out simulations with large negative bias. Use an absolute bias instead, as that metric approaches zero as the simulation gets better.

• The code must return the results of the objective functions as a single tab delimited line. Each objective function result reported in order with no extraneous text.

Step 3: Create an optimization control file

Next you need to create an optimization control file that tells the optimizer and calibrator what it needs to know about your process. You can either start by renaming the Blank_OptiControlFile.txt file that is included with run_scripts in the optimizer repository, or copy one from a previous setup or out of one of the optimizer demonstrations included with the repository.

These instructions assume that you can make use of the AdjustParameters.py code, which works with standard soil and lake parameter files. If you need to work with non-standard file formats or make adjustments to other files, please contact Dr. Cherkauer.

Everything does not need to be defined in the optimization control file yet, as it will be revisited in a later step, instead focus on the follow keywords and sections:

• Parameters

  • Determine how many parameters you are going to include in the optimization process (more parameters will require more simulations to find the Pareto set, while too few will limit the optimizers ability to find a good result).

  • Set N_PARAM to the number of parameters being calibrated

    • N_PARAM <num param> # number of parameters used by the simulation.

  • Add a PARAM keyword for each parameter

    • PARAM <param num> <name> <min val> <max val> <c|m>

    • <param num> indicates which parameter is being defined, from 1 to N_PARAM

    • <name> must be a parameter in the actual file, and must conform to how the soil or lake parameters are named in read_VIC_files.py.

      • To see available paramter names and their spelling use the AdjustParamters.py code with the -h option.

      • AdjustParameters.py foo <Num Layers> <ARC|COL> foo.old foo.new -h <Model Options>

        • <Num Layers> is the number of soil layers in your current application

        • <ARC|COL> indicates the type of soil inputs being used, most people will use text column files (COL), but ARC is available if you are using the ArcInfo input file format

        • Available <Model Options> can be found by typing "AdjustParamters.py" with no command line options. These include flag for Spatial Frost, Groundwater, Lakes, and more.

      • Add "_#" to identify a specific layer (# = 1 to N_LAYERS) parameter to adjust. Note that each layer to be adjusted must be listed separately, the code dose not currently allow for changing multiple layer values the same way.

    • <min val> and <max val> defined the acceptable range for the current parameter

    • <c|m> Use "c" to tell the optimizer that the parameter is to be set to a value in the defined range as a constant. Use "m" to indicate that the existing parameter is to be multiplied by the value in the defined range. Using "c" will set the paramter to the same value for all grid cells, while using "m" will preserve spatial variability by changing the parameter using a multiplier.

• Model Controls

  • Identify the name of the VIC model global file

    • GLOBAL_FILE <file name>

  • Set the name of the model simulation control script

    • RUN_SCRIPT <model run script>

    • Default is set to ControlSimulation.csh, though other names are possible.

  • Define a character string used to separate run time files created by this base simulation, from those created by other simultaneous runs of this program (for example, if you want to calibrate multiple watersheds simultaneously, each should get a separate RUN_ID so that runtime files do not conflict with each other).

    • RUN_ID <model run identifier>

  • Set the name of the return file that MOCOM-UA will seek to obtain the results of each simulation, this file will contain the output from Step 2.

    • RETURN_FILE <file name>

  • Set model options so that the correct inputs and model version are used by the optimizer.

    • MODEL_OPTIONS <options...>

    • This includes things like using a personal copy of the VIC model executable (default is to use the system "vicNl" executable), and identifying major model options (lakes and wetlands, snow bands, groundwater, etc).

    • Available <Model Options> can be found by typing "AdjustParamters.py" with no command line options. These include flag for Spatial Frost, Groundwater, Lakes, and more.

• Objective functions

  • Set the number of tests or objective functions to be assessed, must agree with the output from Step 2.

    • N_TESTS <num tests> # number of test functions to be minimized.

  • Define the name and expected range of the objective functions. These are used to label objective functions in the log file, and can be used when building automated plotting scripts as part of the calibration and optimization process.

    • OBJ_FUNC <test num> <function name> <applied variable> <min val> <max val|inf>

      • <test num> indicates which objective function is being defined, from 1 to N_TEST

      • <function name> is a short name (no spaces) used to identify the objective function

      • <applied variable> is a short name (no spaces) for the variable to which it is applied

      • <min val> and <max val|inf> defined the expected range of objective function values, with "inf" used to indicate that there is no upper limit.

• Queue system requirements

  • Set the name of the SLURM queue to be used

    • SLURM_QUEUE <queue name>

      • Set <queue name> to the appropriate queue for your simulation. For PHIG members, this would normally be set to "phig".

  • Set the length of time required to complete a single VIC model simulation (this is used for the Calibration and Single CPU version of the optimizer, not used for the MPI optimizer).

    • SLURM_WALLTIME <wall time>

      • Format DD-HH:MM:SS, see SLURM manual for more information

• Routing Controls (optional, used by CalculateSimulationStat.py), if you are using the routing model as part of post processing, the following keywords are available:

  • Set the number of routing models that need to be run, typically one for each watershed being routed

    • N_ROUT <num rout>

  • Identify which routing model is going to be run

    • ROUT_TYPE <GIS|ORIG> # GIS for GIS routing model, ORIG for original (Lohman) routing model

  • Setup the command line options for each routing model

    • ROUT_CONTROL <rout num> <rout command line>

    • Command line should include the name (and path if necessary) to the routing model, followed by all command line options.

    • This line will be parsed to identify the names of files required based on the type of routing model being used.

  • Identify the number of routing model output files that will be processed.

    • N_ROUT_OUTPUT <num rout out>

  • Identify the names of each of those output routing files

    • ROUT_OUTPUT <rout out num> <file name>

    • Where <rout out num> identifies the output file for each routing model from 1 to N_ROUT_OUTPUT.

• Comparison files (optional), developed to identify extra files required for comparisons and objective function calculations, these keywords can be used to identify any non-standard file that needs to be copied for each simulation.

  • Set the number of additional files needed

    • N_CMP <num comp files>

  • Set the name of each file

    • CMP_FILE <comp file num> <file name>

Step 4: Setup the simulation control file (RUN_SCRIPT)

Next you need to develop a single overall script that will start the VIC model simulation, wait for it to end and start the post processing, and compute the objective functions when post processing is completed. The script SubmitSimulation.csh, which is part of the run_scripts folder in the repository should have been copied to your base simulation folder in Step 1. This script is already designed to accept the optimization control file and extract from that file information it needs to function. It starts by running the VIC model with the provided global file. It will then run the routing model, if N_ROUT is set. Then it will run the script CalculateSimulationStat.py to compute the objective functions defined in that script.

To customize the script, I suggest that you add your own post processing and objective function calculations between the model simulation and the routing model call, or after the routing model call. Make sure that when complete, your script puts the results of the objective functions into the file name defined in the optimization control file with the keyword RETURN_FILE. That file must contain only the numerical results from the objective functions, in the order defined in the optimization file, with each value on the same line separated by a tab character.

This script must be made executable and must be written in C-Shell.

If you customize this script, it can take only two command line options. First must be the optimization control file, and second must be a flag indicating whether or not the simulation is being run as part of the MPI optimizer (flag = NONE), or is being sent directly to the queue system (flag = queue name). This command line structure is coded into the MOCOM source code, changing it requires change both the single_CPU and MPI code and recompiling.

Test your script by requesting an interactive session, and running all parts of the simulation and analysis with only a single command. Make sure that the name of your script matches that in the optimization control file defined with the RUN_SCRIPT keyword.

Step 5: Test the simulation setup script

Run the command SetupVicSims.sh from within the base simulation folder using the optimization control file you have refined in the previous steps.

This script will create a new file called CopyFiles., where RUN_ID is defined in the optimization control file.

Check this file to make sure that it contains all of the file needed to make your model simulation and post processing work, with the exception of any files defined using absolute paths.

If there are no errors and you are sure that you have defined all of the necessary model option flags in the optimization control file, but still are missing some files then you will need to add those files to the optimization control file using the N_CMP and CMP_FILE keywords. There must be a numbered CMP_FILE line for each file, and N_CMP should equal the total number of files.

When you are done with testing, remove all of the model output and post processing output files to clean the base scenario folder.

Step 6: Test the script base using the calibrator tool

Start by copying the base scenario folder to the Scratch Drive for the cluster system which you will be using for the calibration/optimization process. From this point on, you will start running more simulations and creating a lot more files so you need to the speed and very high storage limits that come with the Scratch Drive system. Anything left on you Home Drive or on the Depot Drive will now become limitations.

The script CalibrateModel.csh reads the optimization control file and then prompts the user to provide values within the range for each of the parameters defined in the control file.

Start the Calibrator and run several simulations while changing the parameter values. You will be prompted after each model simulation to enter new parameters or "Enter" to reuse the previous parameter value (after the first iteration). Use this process to:

• Check that no errors are raised when the submission script is sent to the cluster queue. If errors are raised, review the earlier setup steps. If you followed this process to this point, then new errors are probably:

  1. related to the move to the Scratch Drive so check for broken paths, or

  2. related to the settings for the queue system since this is likely the first time the job has been sent to the queue directly.

• Check that the objective functions are being correctly calculated, stored in the simulation folder under the correct name (defined with RETURN_FILE keyword), and that the values change for each different set of parameters.

If the calibrator passes these tests, then it can be used to facilitate manual calibration of the model. Each new simulation is stored in a new directory (folder name + .####) including the parameters used and the objective function results. It does help to establish a plotting script in the simulation control file so that you can more easily visualize the calibration results.

If the end goal is not a manual calibration, but instead an optimization of the parameters, then move on to the next step in the process.

Step 7: Establish optimization process with the optimization control file

Before starting the optimizer, it is important to adjust some additional parameters within the optimization control file. WARNING: you have two copies of this file now (Scratch Drive and Home/Depot Drive), so make sure that any changes you make to one are also made in the other.

• Define the number of random simulation parameter sets with which to start the test population. The optimizer will generate this number of randomized parameter sets and then select the best test population sized set before starting the optimization process.

  • START <num start | restart file>

    • Set <num start> to a value greater than or equal to the number of sets (N_SETS) in the test population of each generation.

    • If <num start> equals N_SETS than all members are used to start the optimization process.

    • If <num start> is greater than N_SETS then the "best" N_SETS parameter sets are used to jump start the optimizer.

    • If START is set equal to a file named then the optimizer is restarted using the parameter sets and objective function results in the named restart file. Restart files can be generated from an optimization log file (LOG_FILE).

      • Restart files can be used to restart an optimization process that ended early.

      • Addend a '+' to the end of the and the objective function values are not read from the restart file, instead the model simulations are rerun. This allows the optimization of a new model configuration to be started based on the parameter sets chosen by an earlier optimization process.

• Define the size of the test population.

  • N_SETS <num sets>

    • <num sets> is number of simulation parameter sets in the test population. Population size depends a lot on the number of parameters and the number of objective functions used in the calibration. Too few members in the population and the optimizer quickly finds parameter sets that are not dominated by another parameter set. The result is a short optimization process that results in a Pareto set that is not particularly good, since the optimizer did not have to evaluate many generations before it had no further need to look.

    • Too many members in each test population results in long optimization times with many generations and limited advancement late in the process as the optimizer struggles to find a large number of non-dominated parameter sets.

    • Test populations of about 100 members generally work well for ~6 parameters and ~3 objective functions.

    • If using the MPI code, then it is most efficient for N_SETS to be multiples of the number of cores being used. For example, submission to a node with 24 cores should set N_SETS to 24, 48, 72, etc.

• Define a log file which records the steps the optimizer has taken, and the final optimized results.

  • LOG_FILE <optimization log>

• Define variables to control the saving of generational information during the optimization

  • SAVE_SIMS <FALSE|storage directory>

    • If set to FALSE (default) the the optimizer saves only the parameters and objective function results in the LOG_FILE. Actual simulation information will be lost so simulations will have to be redone at the end of the process.

    • If set to a directory then simulation submission scripts and result files will be stored every SAVE_GENERATION.

  • SAVE_GENERATION <num gen>

    • Defines the number of generations between when the simulation submission scripts and simulation results in the SAVE_SIMS directory. Small numbers use a lot of storage space and processor time to copy files, while large numbers mean that potential restart points are further apart.

    • Saving every 25 generations is not a bad compromise.

Step 8: Run the optimization process

There are two versions of the MOCOM-UA optimization program available:

• MOCOM-UA_single_CPU - used for optimizations with a single CPU, or where parallel processing is controlled manually

  • This is an outgrowth of the optimizer developed for the Condor queue system (no longer used on Purdue clusters).

  • MOCOM is run from the command line on a front-end machine, and builds queue submission scripts for each simulation that can be submitted to the queue system individually.

  • The optimization process is primarily completed as a serial process (on a single CPU) where each new simulation cannot be started until the last finishes.

  • By sharing cores within or across nodes, multiple optimization processes can be run simultaneously, but each optimization process will be slower than if run individually on a parallel system.

• MOCOM-UA_MPI - used for optimizations on a parallel computer system using the MPI software package

  • This is the latest version of the code and the one that takes fullest advantage of the cluster computer environment.

  • The MOCOM-UA_MPI code is submitted to a compute node and makes use of all assigned cores to reduce the overall time of a single optimization.

  • To start the optimizer use the following command:

# Here is the general command for submission

sbatch -A [queue] -N (nodes) -n (cores) -t (walltime) [submit file]

# Here is a example of what the submission command looks like

sbatch -A phig -N 1 -n 24 -t 7-00:00:00 Submit_MOCOM-UA_MPI_Job.sh

# where:

# -A phig, requests the 'phig' queue from SLURM

# -N 1, requests a single compute node (1 is typical but can be up to
the maximum nodes available for the queue)

# -n 24, requests 24 cores from the compute node (typically set to the
number of cores per node for the cluster in use)

# -t 7-00:00:00, is the request walltime, typically ask for maximum on
the given queue (this asks for 7 days on SLURM)

# Submit_MOCOM-UA_MPI_Job.sh is a script provided in the repository
that controls the start of the MOCOM-UA_MPI code on the compute node.

• Please check that the variables MOCOM_CORES, MOCOM_EXE and MOCOM_ControlFile are correctly defined at the start of the Submit_MOCOM-UA_MPI_Job.sh script before submitting to the SLURM queue.

Both versions of the model will run independently of the user, so once started you can log out and back into the cluster without affecting the optimization process. Because the single CPU version of the code runs on the front-end system, you will need to record the name of the front end machine (typically --fe01, --fe02, etc.) and log back into that machine to monitor control program activity, submission and completion of individual jobs to the queue can be monitored from any front~end system.

Monitoring the optimizer

The optimization process can take a long time, and it is important to monitor that process to make sure that it is actually moving towards a solution and not simply spinning it wheels. Here are strategies for monitoring the progress of the optimizer:

• Check that simulations are running correctly in the working directory - this is typically named WORK_, where RUN_ID is defined in the optimization log file. This folder exists only when the optimizer is running or if something causes it to fail before the optimization process stops (for example, exceeding the walltime).

  • Folder structure within the WORK directory is as follows:

    • GEN_ggggg

      • GENggggg_mmmm_ss

        • VIC simulation files

  • Where:

    • ggggg is the five digit generation number, starting with 00000 for the random start and 00001 for the first real generation.

      • Old generation folders are removed as they are no longer needed.

      • Members of older generations that are still desirable get copied to the next generation folder before the previous generation folder is removed. These continue to be "parents" of the next generation.

      • Non-desirable members of previous generations are deleted.

    • mmmm is the population member number for the current simulation, from 0000 to Nsets minus the number of members preserved from previous generations.

    • ss is the simplex reflection indicator

      • 00 = normal population member, most common value

      • 01 = member set (mmmm) was better than previous best set, so try additional extrapolation factor of GAMMA

      • 02 = member set (mmmm) is better than the the second largest point, look for an intermediate lower point by doing a one-dimensional contraction

  • Within each folder GENggggg_mmmm_ss is a script called vicsend.csh that when run from within the directory will restart the specific VIC model simulation defined for that folder. DO NOT run this script from within the WORKing folder with the optimizer is running as the folder may cease to exist while your new simulation runs.

• Check that simulations are stored correctly in the save directory - this is defined in the optimization control file using the SAVE_SIMS keyword.

  • If SAVE_SIMS is set to FALSE, no simulation information will be saved outside of the optimization log file, so skip this step.

  • Otherwise, the SAVE_SIMS folder will contain a folder named using the RUN_D keyword, within which will be generational folders for every SAVE_GENERATION, starting with generation

    2. Folder structure is the same as described for the WORKing directory previously.

    3. The vicsend.csh script can be submitted to the SLURM queue to rerun any model simulation.

    4. This folder and its contents will remain even after the optimizer finishes.

• Check the optimization log file - named using the LOG_FILE keyword in the optimization control file.

  • This file contains a running summary of all optimization processes.

  • It will include regular summary tables (based on the SAVE_GENERATION keyword in the optimization control file) that list all parameter sets and objective function results for current generation members.

  • Between those complete summaries, it will tabulate the parameter sets and objective function results for new members of each generation, which can range from a single new member to more that N_SETS new members, if the optimizer decided to run additional set members. Typically, the number of new members will be less than N_SETS for all but the first generation, since some previous generation members will be retained as parents for the latest generation - if the optimization process is successfully moving towards a Pareto set.

  • Check this file early and often to make sure that:

    • parameter selections are staying within your defined ranges,

    • objective function values are being correctly calculated and are changing, and

    • the number of dominated sets is generally decreasing as the optimizer runs.

      • Each parameter set is ranked by objective function performance relative to other parameter sets in the generation. Dominated sets are those that performed less well than other generation members. The optimizer is done when no members of a generation are dominated (this is the Pareto set). In general, the number of dominated sets will decrease as the optimizer progresses. However, this need not be a steady decline. The number of dominated sets can increase substantially when the optimizer finds a new set of parameters with significantly better performance.

      • The optimizer is working if the number of dominated sets generally decreases but gets larger infrequently.

      • The optimizer may be having trouble if the number of dominated sets is staying about the same, generation after generation, as that suggests that it is unable to identify a solution space that is better than the average model performance.

    • If the optimizer appears to be stagnating over many generations, you may need to reevaluate:

      • The selection of parameters - as those being calibrated may not have enough of an effect on the simulation results being evaluated (visible as minimal change in objective function values with large changes in parameter values, so no identified preference for a smaller range of parameters)

      • The range of parameters - as you may be constraining the solution space away from the Pareto set (visible as parameters migration towards and clustering around one edge of their defined range)

      • The sensitivity of the objective functions - as your objective functions may not be providing a sufficient benefit for moving towards a preferred solution (visible as minimal change or random change in objective function values as parameters are changed)

Using the optimizer results

If the optimizer finishes normally, the final table in the optimization log file will be equal to the final or Pareto set. Each of these parameter sets is better than all of the other set members in some way. The final step of optimization process actually requires that you decide, which of the Pareto set members bets suits your needs. This can invovle:

• Selection based only on best performance of a single objective function (may yield poor performance in one or more of the other objective functions)

• Selection based on best overall performance from all metrics (requires compromise, but can yield reasonable results for all metrics)

• Selection based on visualization of the overall results.

I find the best way to deal with the Pareto set is to remove the set members with the most extreme performance criteria (best performance in a single metric, with poor performance in the others) or most extreme parameter values, and then plot simulated versus observed values versus your preferred observed dataset. Then pick the set that looks best to you. This means that the optimization has limited your choices to those it deems best, but you are still using your expert knowledge to select the best of the best, and weed out any members of the Pareto set that are simply unrealistic or do not conform to your notion of what the model results should look like - but that were not well represented by the objective functions selected.

This step is best completed by generating a lot of plots of the results from the Pareto set model simulations.

Application Specific Optimization Files

Here we detail the files and scripts that must be developed for each application. In many cases examples are provided, and may be used to develop similar files for your own application.

Optimization Control File

The optimization control file defines all variables required to run the MOCOM-UA optimizer, and all variables needed by the user defined control script. The file is similar in structure to the VIC model global control file, where there is a key word at the start of each line followed by pertinent information for each key word. A "#" appearing anywhere on a line in the control file indicates that everything after it should be treated as a comment, blank lines in the control file will also be ignored, so feel free to use comments and spaces to clarify the contents of the control file. A demonstration control file is available for download, and is also shown here:

# Define the number of random simulation parameter sets with which to start the
# test population or the name of restart file. If the restart filename is
# followed by a `+`, the model will be reapplied using the restart parameters
# (default uses restart statistics and does not rerun the model

START <num start | restart file>

# Define numbers of each major element

N_SETS <num sets> # number of simulation parameter sets in the test population.
N_PARAM <num param> # number of parameters used by the simulation.
N_TESTS <num tests> # number of test functions to be minimized.
N_CMP <num comp> # number of file names provided for comparisons

# Define the model run script to be used to control the model simulations.

RUN_SCRIPT <model run script>

# Define the VIC model global file

GLOBAL_FILE <VIC global file>

# Define model options that are in use, this includes defining the VIC model
# executable (-V), and options such as the lake model (-L) that require additional
# input files not used by the standard model. Check AdjustParameters.py
# for the latest list of model options.

MODEL_OPTIONS <options list>

# Define a character string used to separate run time files created by this program,
# from those created by other simultaneous runs of this program.

RUN_ID <model run identifier>

# Define a log file which records the steps the optimizer has taken, and the final
# optimized results.

LOG_FILE <optimization log>

# Define variables to control the saving of generational information during the
# optimization, default (SAVE_SIMS = FALSE) is to save only the parameters and
# objective function results in the LOG_FILE. Actual simulation information
# will be lost so simulations will have to be redone at the end of the process.

SAVE_SIMS <FALSE|storage directory> # FALSE or the directory in which to save generational results.
SAVE_GENERATION <num gen> # indicates how often results will be saved.

# Set the name of the return file that MOCOM-UA will seek to obtain the results
# of each simulation.

RETURN_FILE <return file>

# Define all model parameters. There should be one line for each parameter, so
# <param num> should go from 1 to N_PARAM. For each parameter define name, usable
# range and whether it is set as a constant value "c" or through multiplication of
# an existing value "m".

PARAM <param num> <name> <min val> <max val> <c|m>

# Define SLURM queue management.

SLURM_QUEUE <queue name> # name of the SLURM queue to use (e.g., phig)
SLURM_WALLTIME <walltime> # time in DD-HH:MM:SS that a single VIC simulation requires to complete

# Define all routing model requirements. Set the number of routing models 1 to N_ROUT
# that need to be run, typically one for each watershed being routed. These are used
# by the default simulation control script to apply routing models to the current model
# simulation. If not evaluating routed discharge, or writing your own simulation
# control script, these keywords are optional.

N_ROUT <num rout> # number of routing model setups to complete (1 per watershed)
ROUT_TYPE <GIS|ORIG> # GIS for GIS routing model, ORIG for original (Lohman) routing model
ROUT_CONTROL <rout num> <rout command line> # command line (model + options) for each routing model
N_ROUT_OUTPUT <num rout out> # number of output files to expect
ROUT_OUTPUT <rout out num> <file name> # name of each output file

# Define all comparison files. There should be one line for each file, so
# <cmp num> should go from 1 to N_CMP. Define the name of each file containing data
# for the objective function evaluation, e.g. observed discharge. Use these keywords
# to define additional files that must be copied for each model simulation, but that
# are not otherwise identified by the model setup script.

CMP_FILE <cmp num> <obs/true file>

# Define all objective function tests. There should be one line for each test, so
# <test num> should go from 1 to N_TESTS. For each test define objective function
# name (e.g., variable function is applied to and range of resulting values (use
# "inf" to indicate no maximum limit).

OBJ_FUNC <test num> <function name> <applied variable> <min val> <max val|inf>

# Define all additional lines, these will not be used by MOCOM-UA but will instead
# be passed directly onto the simulation control script, thus they should be defined
# to fulfill any functionality needed by the control script. Do not use any of the
# key words already defined. There is no limit to the number of lines that follow.

ControlSimulation.csh

• This script should be in the LOCAL optimization directory.

• This program should be edited to call the user developed programs for computing objective functions and plotting the current results. Plotting is optional.

Process relies on the following additional libraries and scripts installed on the brown cluster as part of the /depot/phig/apps library:

• AdjustParameters.py - adjusts soil parameters for the simulation set.

• read_VIC_files.py - reads and writes a number of standard VIC model input and output file formats.

If running on a system other than the brown cluster, these libraries will need to be installed and included in your PYTHONPATH variable.

submit.job

• This script should be in the LOCAL optimization directory.

• Wrapper that submits the entire optimization process to the PBS queue.

Application Specific Analysis Script

• This should be an executable program in a local path that is included in your PATH environment variable (e.g., ~/bin)

• This is a custom file created by the user that computes the values of the selected objective functions once the model simulation is complete.

• Should return objective function results to stdout, on a single line, in the order in which they appear in the optimization control file.

• Objective functions will be minimized by MOCOM-UA, so best value must be the smallest value.

  • Some objective functions must be inverted for this to work.

  • For example, Nash-Sutcliffe Efficiency (NSE) provides values between negative infinity and 1, where 1 is a perfect simulation. To use NSE with the optimization, the analysis script should return the value "1 minus NSE" as that will go to 0 from positive infinity as the simulation improves.

• It can be designed to make use of any information in the Simulation Control File.

Need to develop a sample program that reads more of the information out of the control file, and therefore requires fewer modifications.

Application Specific Plotting Script (Optional)

• This should be an executable program in a local path that is included in your PATH environment variable (e.g., ~/bin)

• This is a custom file generated by the user to plot simulation output, typically versus observed values.

• Output from this script serves as a way to visually inspected simulation data, both when using CalibrateModel.csh and during the optimization process.

• It is not a required part of the optimization process, but as the generated plots can be saved as part of the simulation directory, they can be viewed at anytime during or after they are generated for a visual check on optimization performance.

Simulation Run Script

The Simulation Run Script (RUN_SCRIPT in the optimization control file) is the key to making the optimizer work. Each simulation run script must:

• Be an executable script or program that MOCOM-UA can execute using a simple system call.

• The simulation run file must accept one command line argument, the name of the Simulation Control File (see description under General Programs).

• It must be placed in a location include in your PATH environment variable.

If you are optimizing a VIC model simulation, then you can probably make use of the RUN_SCRIPT, SubmitSimulation.csh.

• It should be usable for most VIC model optimization applications, without modification.

• Everything it needs to control the simulation should be defined in the Optimization Control File.

• This is stored in a common directory and should not need to be edited.

• This script builds a local version of "vicsend.csh", which wraps the submission of the ControlSimulation.csh script appropriately for submission to the PBS queue.

Note

To Fix

The application specific analysis script, the application specific plotting script and the the simulation walltime are all hardwired into this version of the script, so the script must be updated to use the variables defined in the Optimization Control File.

If you are trying to optimize something other than the VIC model, then you should write your own Simulation Run Script.

• The script is called by MOCOM-UA and provided a version of the Simulation Control File that defines the parameter settings for the current simulation, and passes along any additional information from the Optimization Control File that was not used by MOCOM-UA.

• Typically this information is put at the end of the Optimization Control File, but it can appear anywhere in the file as long as it does not reuse keywords required by MOCOM-UA.

• Test the script using a sample Simulation Control File before starting the optimizer.

General Optimization Files

Here are details related to the files and scripts that must be installed or otherwise available before running an optimization. For most applications, these will not need to be modified by the user.

MOCOM-UA program

The multi-objective optimization scheme used by this program is the Multi-Objective Complex Evolution - University of Arizona (MOCOM-UA) method (Yapo, et al., 1998). This is based on the Shuffled Complex Evolution - UA (SCE-UA) single-objective method but was modified to work with multiple objective functions. This program was written by Keith Cherkauer using the description of the algorithm provided by Yapo, et al. (1998) and does not make use of any source code from UA. It was specifically designed to be very general so that it could work with any script that given a control file with some basic parameters and the potential for unlimited user added information, will return the results of multiple objective functions. It should be noted that all objective functions must result in lower values as the simulated and comparison files get closer together.

General description of MOCOM-UA optimization process

• The MOCOM-UA program is initialized by generating START parameter sets randomly based on the defined range of parameters in the control file. These parameter sets are used sent to RUN_SCRIPT and objective functions are evaluated on the simulation results and returned to MOCOM-UA in RETURN_FILE.

• Next the simulations are ranked starting by finding the globally non-dominated sets (rank 1). Members of rank 2 are dominated only by members of Rank 1, this continues until all memebers of the current population have been ranked. For a discussion of non-dominated sets, please see Goldberg, "Genetic Algorithms in Search, Optimization, and Machine LEarning" (1989).

• Sets in the current generation are then sorted by rank and the optimizer starts working with the N_SET highest ranked members of the starting generation. If N_SET is less than START, the optimizer will start using the N_SET lowest ranked parameter sets from the START members of the initial generation.

• Next the probability of each population member reproducing is calculated where the lowest ranks have the highest probability of reproducing.

• Now each member of the highest rank is bred with N_PARAM randomly selected parameter sets from the lower ranking population using their probability of reproduction.

• The N_PARAM+1 members of the breeding population are fed to amoeba, which is a simplex optimization scheme from Numerical Recipes. The amoeba algorithm will try reflection and contraction using one or more of the lower ranking members and the highest ranked member in order to find a new parameter set that results in improved objective function solutions. If it fails, MOCOM-UA tries again by selecting a new breeding population (Step 5), until a better set of parameters has been found or the maximum number of iterations has been exceeded.

• If successful, each member of the highest rank in the previous generation has now been replaced by offspring combining features of it and N_PARAM better parameter sets.

• Rank is now computed for all members of this new generation, generation members are sorted by rank and then new reproduction probabilities are calculated.

• If the number of ranks still exceeds the maximum number required to end the optimization (Rmax = 1 for a complete optimization), then the optimizer returns to step 5.

• If the number of ranks in the current generation are less than or equal to Rmax, then the optimizer will output the final results of the optimization. If the number of Ranks in the final generation equals 1, then the members of the generation define the Pareto set, i.e. none of the members are dominated by any other so each set of parameters is just as likely as all of the others to be the optimal solution of the simulation model.

Availability:

• On conte and other RCAC cluster systems, the latest versions of this program are compiled and available to run in the shared directory: /depot/phig/apps/VIC, just make sure that this directory appears in your PATH environment variable.

• The source code can be accessed here for installation on other systems.

Versions:

• MOCOM-UA_single_CPU.c - designed for use on a single CPU. This is the default version of the optimization code, and should work on all platforms on which the VIC model can be run.

• MOCOM-UA_MPI.c - Still under development. Should allow the optimization routine to take full advantage of parallel processors to more rapidly complete the calibration process.

• MOCOM-UA_newer.c - Depreciated. Original version of the code, has been replaced with MOCOM-UA_single_CPU.c but still available for backwards comparability with setup files.

Dependencies:

• None.

Simulation Control File

The Simulation Control File is produced by MOCOM-UA using the contents of the Optimization Control File to define parameters specific to a single iteration of the Simulation Run Script. There should be no reason for the user to modify the structure of this file as tit would require editing and recompiling the MOCOM-UA program. Note that information at the end of the Optimization Control File (and not used by MOCOM-UA) is passed directly to this file. Therefore, it is available for use by both the Application Analysis Script and the Application Plotting Script (defined under Application Specific Optimization Files).

The following is an example of a simulation control file, note that most of the comments embedded in the Optimization Control File will be lost as the file is constructed by MOCOM-UA. Some of the differences between this and the Optimization Control File include: The definition of OUT_DIR the directory created to hold the simulation and the location of the Simulation Control File. The list of parameters now contains specific parameter values rather than a range, these values must be used by the Simulation Run Script to modify the performance of the simulation resulting in changes to the objective functions. After the list of comparison files, most of the demonstration file contents were copied directly from the Optimization Control File, so the purpose of those later lines is simply to provide required information to your Simulation Run Script that is not not significant to MOCOM-UA but is used by the script.

OUT_DIR WORK_VicModel_Try4/GEN00038//GEN00038_0000_00/
RUN_ID VicModel_Try4
N_PARAM 9
PARAM 1 Ds 0.065013 c
PARAM 2 Dsmax 0.973148 c
PARAM 3 b_infilt 0.016001 c
PARAM 4 expt_1 19.139690 c
PARAM 5 expt_2 17.752407 c
PARAM 6 expt_3 9.844997 c
PARAM 7 bubble_1 25.549908 c
PARAM 8 bubble_2 43.885036 c
PARAM 9 bubble_3 17.240543 c
RETURN_FILE VicModel_Try4.results
N_CMP 2
CMP_FILE 1 SEPAC_Plots_drainflow_Drainflow.86_08.20m.NoData.dly
CMP_FILE 2 SEPAC_Plots_wtable_WTavg_1988_90_dly.txt
N_SIM 2 # number of routing control files
CMP_COL 1 3 # column of subsurface drainage flows
CMP_COL 2 4 # column of water table depths
CMP_MULT 1 10.0 # multiplier for subsurface drainage flows
CMP_MULT 2 1.0 # multiplier for water table depths
SIM_MULT 1 1.0 # multiplier for subsurface drainage flows
SIM_MULT 2 100.0 # multiplier for water table depths
GLOBAL_FILE SEPAC_InputFiles/global_4.1.2_DWM_NewForcing.opti # global file for VIC model simulation
SIM_OUTPUT 1 drainage_39.0625_-85.5625 4 1.0 # name of routing model output file
SIM_OUTPUT 2 drainage_39.0625_-85.5625 7 1.0 # name of routing model output file
SIM_COL 1 4 # column of drain flow
SIM_COL 2 7 # column of water table depth
SLURM_QUEUE phig
SLURM_WALLTIME 00:30:00

AdjustParameters.py

This Python script reads VIC model soil files and adjusts parameters according to the values in the Simulation Control File generated by MOCOM-UA as part of the optimization process. It is also used by CalibrateModel.csh to adjust parameters for each calibration run.

Availability:

• On conte and other RCAC cluster systems, the latest versions of this program are compiled and available to run in the shared directory: /depot/phig/apps/script/VIC/, just make sure that this directory appears in your PATH environment variable.

• The source code can be accessed here for installation on other systems.

Versions:

• As of 18 Aug 2016 only one version is available.

Dependencies:

• This program requires the Python libraries:

  • read_VIC_files.py

Other Files

This section discusses related, but not required file and scripts that can be used as part of the calibration of the model.

CalibrateModel.csh

• Runs the SubmitModel.csh without entering MOCOM-UA, so runs the model once with a selected set of parameters.

• This program should not need to be modified by the user, all it needs to know about the current job is in the control file.

• It is a useful tool to test that your model simulations are running and that objective funtions and plotting scripts are working correctly before starting the optimization, since this program uses all of the same programs and files.

• This assumes that there is a plotting script generating a postscript plot file in the current directory. Might want to allow it more flexability in the future.

Python Library Modules

There are a number of python library modules that must be installed in your PYTHONPATH in order to make the optimization programs described above work. These libraries include:

• read_VIC_files.py (/depot/phig/apps/source/python/lib/read_VIC_files.py); there is a copy in /depot/phig/data/Projects/OPTIMIZATION/read_VIC_files.py that should be deleted

• datefuncs.py (/depot/phig/apps/source/python/lib/datefuncs.py); there is a slightly older version in /depot/phigs/apps/script/lib/datefuncs/py that should be deleted. This should be depreciated in favor of using Python datetime libraries.

• read_arcinfo_files.py (/depot/phig/apps/source/python/lib/read_arcinfo_files.py)

• stats.py (/depot/phig/apps/source/python/lib/stats.py); there is a newer copy in /depot/phig/data/Projects/OPTMIZATION/stats.py - not sure what was changed.

How to optimize a model simulation

The following instructions are for optimizing the VIC model on an RCAC maintained cluster computer.

• Running on different systems requires having the optimization code installed (see above for General Optimization Files and Other Files for information on what files and programs must installed locally).

• To optimizing a different model, you will need to edit SubmitSimulation.csh

• Start by making sure that you have all of the input files required to run your model simulation.

• Next create the Application Specific Analysis Script that will calculate your objective functions for each simulation.

  1. This script will need to know what the model output file(s) are and where the observed calibration data file(s) are.

  2. It will need to be able to run from anywhere when you are logged into the system.

  3. Your objective functions MUST be minimized as the simulation gets closer to the observations.

  4. It will need to return the objective function values as a list of values in the same order as they are defined in the Optimization Control File.

• If you want to view your results during and after the optimization process, without running a new plotting script, you should create an Application Specific Plotting Script.

  1. This script will need to know what the model output file(s) are and where the observed calibration data file(s) are.

  2. It will need to be able to run from anywhere when you are logged into the system.

  3. You can use whatever plotting program you wnat, as long as it can be run from the command line or from within the RUN_SCRIPT.

     1. I suggest NOT using MATLAB or other plotting programs with a large overhead, as you do not want to significantly impact your simulation times by repeatedly loading large plotting software packages into memory.

     2. Both GMT and the Python library plotmatlib are fairly efficient for generating figures quickly.

  4. Currently the CalibrateModel.csh script looks for an EPS file in the simulation directory, and will display that between calibration sets.

• Create the Optimization Control File.

  1. This will define the scale of the optimization and the number and range of the parameters used.

  2. It should also direct the RUN_SCRIPT to the Application Specific Analysis Script and the Application Specific Plotting Script.

  3. Optimization parameters to pay attention to:

     1. START: As this will be the first time you run an optimization, you should set START to a large value for the random start process, typically several times the value of N_SETS

     2. N_SETS: Number of parameter sets to maintain in the current generation. This MUST be larger than the number of parameters.

The output file <submit script>.o<jobID> will generally include the warning:

Warning: no access to tty (Bad file descriptor).

Thus no job control in this shell.

This warning reminds you that this is a batch job that cannot be acted upon. If you see this warning, then your job was submitted and at least started correctly on the PBS queue.

Note

From Laura:

Files needed for optimization. We are trying to balance the need to maintain code provenance, optimizer function (some code needs to be available no matter what generation sub-directory the optimizer is running in and efficient read/write times.

All code in the central repository (/depot/phig/...) should remain there. That is will the centralized code will be maintained. Some needs to be copied locally, as follows:

a) PYTHON scripts that are also used outside of optimization should remain central (maybe all moved to : /depot/phig/apps/source/python/lib/ ??) AdjustParameters.py (/depot/phig/apps/script/VIC/AdjustParameters.py)

read_VIC_files.py (/depot/phig/apps/source/python/lib/read_VIC_files.py; there is a copy in /depot/phig/data/Projects/OPTIMIZATION/read_VIC_files.py that should be deleted)) datefuncs.py (/depot/phig/apps/source/python/lib/datefuncs.py; there is a slightly older version in /depot/phigs/apps/script/lib/datefuncs/py that should be deleted) read_arcinfo_files.py (cannot find on /depot/phig) stats.py (apps/source/puthon/lib/stats.py although there is a newer copy in /depot/phig/data/Projects/OPTMIZATION/stats.py)

b) Official optimization run scripts, that could be run from the central repository (maybe move all to /depot/phig/apps/VIC/bin/ ?)

• SubmitSimulation.csh (/depot/phig/apps/script/VIC/SubmitSimulation.csh)

• MOCUM-UA_single_CPU source code and executable (/depot/phig/apps/VIC/MOCOM-UA_single_CPU)

• CalibrateModel.csh (not currently available on /depot/phig got copy from the optimization tar file)

c) Optimization run scripts, that should be copied from the repository because they have to be customized (maybe move all to /depot/phig/apps/VIC/ ?)

• ControlSimulation.csh (future version can be automated so editing is not needed) (no communal copy, but three copies exist in /depot/phig/data/Projects/OPTIMIZATION/)

• submit.job (this is a stupid name, needs to change; no centralized version exists, only in /depot/phig/data/Projects/OPTIMIZATION/)

• SetupVICOptimization.csh (future version can be automated so editing is not needed) (haven't found on /depot/phig/ probably have a typo) copyfiles.csh (/depot/phig/apps/script/VIC/copyfiles.csh; another copy exists in /data/Projects/OPTIMIZATION/)

d) Sample files that are available, but need to be copied to local directory and customized

• MOCUM-UA control file - An example, e.g. SEPAC_try2_MOCUM-UA_control_file.asc, is in the tar file Current VIC global file, veg param, veg library, soil param and forcing files.

• Observed drainflow file

• Program to calculate objective functions.

• Plotting program to make graph of simulated versus observed drainflow to be saved in the output (optional).

• VIC executable

Note

Naming Conventions

Each simulation is stored in a self-contained directory that includes everything required to rerun it exactly as done by the optimization program.

During the optimization process:

• Simulations run for the current generation will be stored in a temporary working directory: WORK_<model run identifier>, where <model run identifier> is defined in the control file.

• Within the temporary working directory will be another subdirectory named: GEN<generation number>, where <generation number> identifies the current optimization generation. The starting generation is Number 0.

• Within the generation specific folder will be additional subdirectories named: GEN<generation number>_<setcnt>_<type>, where

  • <generation number> should be the same as the parent directory,

  • <setcnt> is the ID number of the current member of the generation (starts with 0, and should reach NUM_SETS when generation is completed), and

  • <type> indicates, which step of the optimization process was reached. The following values are possible based on what the amoeba function tries to do to improve the results:

    • 0 = reflected simplex from the highest point, this always happens other types are only used if this result is better than existing best parameter set

    • 1 = extrapolate past simplex reflection by GAMMA

    • 2(+i) = solution from try 0 is better than 2nd best point, try contracting the reflected point, each point tried increases value of i.

  • For each generation, the optimizer only generates new members from the worst (highest ranked) members of the current solution set.

    • It selects each member of the highest rank set, and NUM_PARAM additional sets of parameters selected randomly from the current generation and the amoeba function works with those to select a new set of parameters that will provide a result better than that from the parameters from the righest rank member.

    • For example, if there are 3 members of the worst rank in Generation 1, then Generation 1 will only generate 3 new members, with the following exceptions:

      • If the amoeba function decides that the original parameter set can be improved on, then there will be additional folders for each attempt amoeba makes to improve on the solution.

      • If the amoeba function fails to find a better solution (including if parameters are found to be out of bounds or the model does not solve after several tries), then it will try the entire amoeba process again. So there will be additional folders.

    • Therefore the number of files in each generation directory in the working folder can have a different number of sub-directories.

Saved optimization generation files:

• Simulations representing the best solutions of the objective functions are stored in the directory SAVE_SIMS at an interval defined by SAVE_GENERATIONS. Both parameters are defined in the control file.

• The contents of the temporary working directory for the current generation is not the same thing as the save generation.

  • The temporary working directory can have multiple members of the new generation that are not improvements over members of previous generations.

  • Only the best parameter sets are preserved for the next generation.

  • This means that members of older generations can continue to exist in the current save generation because they represent objective function results that no later parameter set has been able to improve.

  • Because of this, saving any particular generation (e.g. that defined by SAVE_GENERATION, or the final generation) requires a more complex process than a straight copy of the current working directory.

• The process by which the current save generation is preserved is described here:

  • Within the directory SAVE_SIMS, the optimizer will open a new directory named for the next save generation and copy the contents of the last save directory into it. For example, the random start generation is saved to "GEN_00000", and if SAVE_GENERATIONS is set to 10, then immediately after saving the random start generation the optimizer will open a new directory called "GEN_00010", and copy all sub-directories from "GEN_00000" into the new directory.

  • The optimizer will then run all simulations to evaluate the first generation.

  • When it is done, it will identify which members of the first generation are improvements over members of the starting generation. These sub-directories will be copied to "GEN_00010", while the members of the starting generation that have been replaced will have their sub directories removed.

  • At the end of the first generation, the "GEN_00010" directory down SAVE_SIMS will contain a mixture of files named "GEN00000_..." and "GEN00010_...". There will only be NUM_SETS sub-directories in the save directory.

  • When the the optimizer has completed the tenth generation, the directory /GEN_00010 will contain a mixture of sub-directory names that represent the parameter sets that make up the best simulation possibilities after ten iterations of the optimizer. These will be copied to the new directory GEN_00020, and the contents of GEN_00010 will not be changed further, thus preserving the results, plots, analysis and files required to recreate that entire generation.

Restart Process

The optimization process can be restarted from any saved generation. This allows you to recover from a failure of the original optimization (system crash, exceeding walltime, etc.), while preserving the status of the optimization process. It also allows you to run a staged optimization, where a first optimization uses a simplified model or coarser time step, then the restarted optimization can use a full model simulation to refine results.

Note

Not sure that this latter processes is effective, but it was used for calibrating in energy balance mode then adding soil frost later. The hypothesis was that a minor tweak to the model should leave the original solution space in a similar region. I will have to see if I can dig up anything on those original test, but I think that the model solution space is not nearly that consistent.

To restart an optimization process so that it picks up where it left off (so no change in model configuration), the following steps should be completed in the same directory used to conduct the original optimization:

• Copy the LOG_FILE (defined in the optimization control file) into a new file named to indicate that it is for restarting the optimization run.

• Open the new file in the editor.

• Search for the last full saved generation - searching for the text "Current generation" should skip over the intermediate unsaved generations.

• Saved generations should look like the following in the log file:

  Current generation for try 30:
  ==========
  Ds Dsmax b_infilt expt_1 expt_2 expt_3 bubble_1 bubble_2 bubble_3 tanB = ( IOA fract PBI fract MAE fract ) Rank Name
  0: 0.016733 1.7112 0.34544 18.754 10.206 8.8939 56.819 97.046 36.116 0.030229 = ( 0.34359 1.9054 0.26695 ) 1 GEN00028_0000_00
  1: 0.015336 2.0021 0.38943 16.397 9.23 10.35 64.976 80.478 53.17 0.036705 = ( 0.34396 0.001457 0.000204 ) 1 GEN00027_0000_01
  2: 0.0058379 0.98813 0.38905 16.854 10.65 12.748 60.35 57.492 64.259 0.032101 = ( 0.34448 4.4491 0.62333 ) 2 GEN00028_0003_00
  3: 0.029781 1.1746 0.3327 14.7 12.723 13.911 61.17 66.734 63.393 0.040137 = ( 0.35711 0.21061 0.029508 ) 2 GEN00013_0000_12
  4: 0.025459 1.9441 0.33334 15.346 11.764 12.849 63.814 75.19 52.609 0.039043 = ( 0.35727 0.10719 0.015017 ) 2 GEN00029_0002_02
  5: 0.055098 2.4755 0.26345 16.626 12.648 13.717 53.219 82.689 38.955 0.041233 = ( 0.35414 0.6739 0.094416 ) 2 GEN00018_0003_02
  ...
  
  <NUM_SETS>: ....
  ==========
  

• The header lists the parameters in the order that they are defined in the optimization control file, followed by the objective functions in "()".

• Each row until the final "=========" provides the parameters and objective function results for each simulation within the current generation. There should be NUM_SETS rows under the header.

• To convert this optimization file into a restart file:

  • Delete all lines in the restart file outside of the records for the last saved generation as marked above.

    • The first in the file must be "Current generation...", followed by three additional lines of header material.

    • The resulting file must be at least NUM_SETS + 4 header lines long. Additional lines past the end will be ignored.

• Next copy the optimization control file to a new name that indicates it is being used for a restart of the optimization.

• Edit the entry after the keyword "START".

  • To restart the optimizer, change the number representing the number of random start simulations to run, to the file name of the restart file you just finished creating.

  • If you have made changes to the model or it options, which will render the original objective function values inaccurate, then add a "+" to the end of the restart file name. This will force the optimizer to rerun all of the members of the restart generation so that it updates all objective function values before continuing the optimization process.

• If you were using the "submit.job" script, and created a new restart file, then you will need to edit the submit.job script to use the restart control file and not the original control file.

• You can now use the restart control file to jump start the optimizer from the last generation.

To restart the optimization using a new model configuration, such that objective function results are not expected to match the original. You should follow the above instructions in a new directory, or change the RUN_ID, and optimization log file name and run in the same directory. If you do not then restarting the optimizer will overwrite your original optimization results.

Development Notes

• Have generated a methodology for passing a control file between the MOCOM-UA program the the function that will be evaluating whatever the models is. It will need to include information about

  1. parameter values,

  2. generation number,

  3. whether or not data is to be saved after simulation,

  4. name of the expected objective function evaluation file,

  5. name of the comparison file,

  6. perhaps information concerning the run location of the main program and storage location for saved data?

• Need to finish the single CPU version first, before starting to dig into the MPI programming required to get MOCOM-UA running on the cluster systems. For MPI programming, amoeba will need to be split across different CPUs. Amoeba may have to have to run the model several times, but each call to amoeba is for a specific member of the current generation's population.

• Current MOCOM-UA version seems to have collected a couple of bad parameter sets (all parameters are 0) and seems unable to discard them. Need to update code to look for failed parameter sets and replace them with valid values (either from previous generations or simply by generating them randomly). This is becoming a larger problem!!!!

• Might want to add information about the objective functions to the setup control file. Information such as function names and ranges would be helpful in generating log files and plotting optimization results. Including expected minimum and maximum values would also be useful.

• I call too many things "try" in MOCOM-UA, need to develop a clearer vocabulary.

Development To Do List:

• Dr. Cherkauer will clean up directory structure for source code

• Update MOCOM-UA to make use of multiple CPUs

• Complete documentation for basic optimization on Conte

• Bundle complete optimization example and link to this page.

Notes Related to the autocal_only version of the model

In this version of autocal, nothing is kept in the autocal code directory except the codes and copies of sample run scripts needed to run the approach.

The run_calib/ directory is copied to somewhere where it's easy to navigate between the calibration effort and the model's regular runs. It contains scripts that run MOCOM and the local model implemementation.

The latest modifications to the code were to take PCIC's upgraded version and add back some of the functionality to store/keep intermediate run outputs (rather than just output a log file listing final parameter values). Mainly for convenience -- with the run outputs available, any run from the process could be easily plotted/analyzed, and the resulting parameter files copied directly into the model parameter directories for use if selected. For troublesome some calibrations, it also allowed for more diagnosis of the progress of the optimization.

AWW-Jan2016

mocom-ua.aww was the original code from UW as modified & upgraded by A. Wood. PCIC later went through the code and report to have fixed several bugs, though this is unverified by code review. The arguments to the code change differ slightly, and PCIC added a daily stats calculator so that if daily flows are available, they could also be used.

=================

Two more recent versions of MOCOM sent by PCIC/Katrina Bennett around 2009. The apr-dm one seems most recent.

PCIC apparently found some bug fixes that were important, so this should be used.

The src_bin/ dir has the apr-dm/ codes

=======================

the PCIC version removed some functionality which allowed the user to store intermediate run results, so

Optimization Test Sets

The following tar files contain control files, data sets and scripts that have been used for the optimization of different scenarios, mostly of the VIC model, though other cases can also be provided. When posting an optimization test package please make sure that all files required to run the optimization are provided.

• Sinusoidal evaluation - This package contains a simple code to produce sinusoidal output given two amplitudes and two frequencies. It was designed to provide a simple and quick test case for the optimization code.

• Streamflow adjustment - This package contains a program that will adjust observed streamflow using ### parameters. It was designed to provide a more realistic test for the optimizer than the sinusoidal function generator.

• Subsurface drainage evaluation - This package contains data files for the SEPAC subsurface tile drained fields and was used to identify optimized parameters for Sarah Rutkowski's Masters research. The SEPAC site was run as a point scale using the modified VIC model with water table, subsurface drainage and urban algorithms. Output subsurface drain flow and water table depth were evaluated versus observations for a 20 m drain spacing using Nash-Sutcliffe Efficiencies (inverted, e.g. 1. - NSE) and Mean Absolute Error (MAE) objective functions for both. Drainage was excluded from the calculation when both simulated and observed drainage was zero, while water table depth was excluded when the observed depth was greater than 120 cm (such measurements were considered suspect by Dr. Kladivko who oversaw data collection).

Package Contents

Download the package from the hydrodat ftp site using this link.

ControlSimulation.csh (Jan, 2011) - Simple script designed to control execution of the VIC model for a given global file. Does not currently rout runoff and baseflow to streamflow, and was not designed to change model parameters.

MOCOM-UA_newer.c (Jan, 2011) - This program uses the MOCOM-UA multi-objective optimizing scheme to optimize anything included in the optimization script that returns statistics that can be minimized to find the 'actual' solution. The result will be a set of simulations that define the Pareto solution set. Given a sufficiently large population, the Pareto set should define the best set of solutions obtainable with the calibration script.

RelateNijssenAndStandardVicBaseflowParams.py - This is not really a program, it appears to just be the equations for relating Bart Nijssen's baseflow parameters to the more traditional ones used by the VIC model. Both are acceptable as inputs to the VIC model, the Nijssen parameters remove the relationship between baseflow curve shape and the depth of the bottom soil layer, which should make it a more robust set of calibration parameters. Baseflow is still dependent on the depth of the bottom layer, but with the Nijssen parameters that relationship is computed within the VIC model so changing the bottom layer depth will not affect the shape calibrated with the given baseflow paramters.

SubmitSimulation.csh (Jan, 2011) - This script is designed to take a list of parameter adjustments from MOCOM-UA and farm out VIC model simulations for each. It will then wait for all simulations to finish, gather simulation performance statistics and return those to MOCOM-UA. At least that is the hope, have to determine how far I got with this script before knowing how useful it is. Also need to check radon source code, since I believe that I have a working version of the basic submission script.

check_optimization.scr (May, 2003) - This is an example script for running the VIC model within the MOCOM-UA_newer wrapper, it has not been modified to work with multiple CPUs but should still work for a single CPU though I do not believe that it has been run since moving to Purdue. This script was written to run all processing steps necessary to create new parameter files, run the model, and process the output for a new calibration parameter set.

check_params.txt (April, 2003) - Sample parameter definition file. Establishes acceptable ranges for parameters being used by MOCOM-UA. Currently MOCOM-UA uses values drawn randomly from a linear fit between the values, but it might be valuable to allow it to use a logarithmic or power function, since some parameters range over several orders of magnitude.

compute_daily_weekly_monthly_runoff_stats.c (May, 2003) - Modified from compute_new_weekly_runoff_stats.c, this program computes the mean squared error, peak difference and negative number of sign changes between the simulated and observed discharge records. Discharge is converted from daily to weekly and monthly values and test statistics are computed for all time steps beginning with the given starting date. Not tested since moving to Purdue, but no real reason that it would not work.

create_restart_file.py (May, 2003) - This script reads an optimization log produced by MOCOM-UA_newer and strips out the parameters, statistics, rank and solution set number for the selected generation of the optimization record. The resulting restart file can be used to restart MOCOM-UA_newer at the given generation. This is useful for restarting an optimization run that was killed prematurely as well as to start a new optimization using a given set of parameters (e.g. using VIC in WB mode to obtain a Pareto set and then using VIC in EWB-FS mode to refine that set). Not tested since moving to Purdue, but no real reason that it would not work.

• Modified to remove reading and writing of simulation number since it is no longer included in the log file.

• Need to update this program to work with optimization control file. Should be able to get everything it needs from the file, and to create a new control file for the optimization process.

modify_flow.py (May, 2003) - Not sure what I was doing with this, though it might have been to modify an observed time series in a controlled fashion in order to test the MOCOM-UA program, since in this case the optimal solution should be to minimize the changes to the observed time series.

Optimization Demonstration

Based on optimizations conducted by Sara Rutkowski as part of her Master's thesis work. The optimization focuses on simulated subsurface drainage versus observations of subsurface drainage at the SEPAC field site. This tutorial assume that you have a version of the VIC model source code installed and compiled that include subsurface drainage.

1. Download the file SubsurfaceDrainageOptimization.tar and unpack into a working directory.

   1. This is now included in the GitHUB resposity, at https://github.rcac.purdue.edu/PHIGOrganization/Optimization/tree/master/Purdue/TestDatasets

2. Download and compile MOCOM-UA_single_CPU.c for your current platform.

References

Yapo, P., H. Gupta, and S. Sorooshian, 1998: Multi-objective global optimization for hydrologic models, J. Hydrology, 204(1), 83-97.