SWAT+ Calibration Example using pySWATPlus and pymoo¶

This notebook demonstrates how to use the pySWATPlus package to calibrate SWAT+ model parameters using the pymoo optimization library. The example focuses on calibrating parameters in the plants.plt file and minimizing a user-defined objective function.

---

1. Import Required Libraries¶

First, we import the necessary libraries, including pySWATPlus for SWAT+ model interaction and pymoo for optimization.

In [16]:

import numpy as np
from pySWATPlus.TxtinoutReader import TxtinoutReader
from pySWATPlus import SWATProblem, minimize_pymoo
from pymoo.algorithms.soo.nonconvex.cmaes import CMAES
from pymoo.termination import get_termination
from pymoo.util.normalization import denormalize

import numpy as np from pySWATPlus.TxtinoutReader import TxtinoutReader from pySWATPlus import SWATProblem, minimize_pymoo from pymoo.algorithms.soo.nonconvex.cmaes import CMAES from pymoo.termination import get_termination from pymoo.util.normalization import denormalize

---

2. Define the Objective Function¶

The objective function is a user-defined function that evaluates the performance of the SWAT+ model for a given set of parameters. It must:

• Accept a single dictionary argument containing the calibration parameters.
• Run the SWAT+ model with the provided parameters.
• Calculate and return an error metric based on the model's output.

In [17]:

def function_to_minimize(dict_of_params):
    """
    Objective function to minimize. It runs the SWAT+ model with the provided parameters and returns an error metric.

    Parameters:
        dict_of_params (dict): A dictionary containing the calibration parameters and other necessary information.
                               Must include the key 'calibration_params' with the format:
                               {filename: (id_col, [(id, col, value)])}

    Returns:
        Tuple[int, Dict[str, str]]: A tuple containing the error metric and a dictionary with the simulation results.
    """
    # Extract calibration parameters and path to the SWAT+ TxtInOut folder
    calibration_params = dict_of_params['calibration_params']
    path_to_txtinout = dict_of_params['path_to_txtinout']

    # Initialize the TxtinoutReader and copy the SWAT+ project to a temporary directory
    reader = TxtinoutReader(path_to_txtinout)
    tmp_path = reader.copy_swat(dir=None)  # Copy to a temporary directory
    reader = TxtinoutReader(tmp_path)

    # Run SWAT+ with the provided calibration parameters
    txt_in_out_result = reader.run_swat(calibration_params, show_output=False)

    # Initialize a new TxtinoutReader to read the results
    result_reader = TxtinoutReader(txt_in_out_result)

    """
    The following steps should include:
    1. Reading the simulation results.
    2. Gathering observed data.
    3. Calculating the error metric based on the observed and simulated data.
    """

    # For demonstration, we return a random error metric
    rng = np.random.default_rng()
    return (rng.random(), {'test_calibration': result_reader.root_folder})

def function_to_minimize(dict_of_params): """ Objective function to minimize. It runs the SWAT+ model with the provided parameters and returns an error metric. Parameters: dict_of_params (dict): A dictionary containing the calibration parameters and other necessary information. Must include the key 'calibration_params' with the format: {filename: (id_col, [(id, col, value)])} Returns: Tuple[int, Dict[str, str]]: A tuple containing the error metric and a dictionary with the simulation results. """ # Extract calibration parameters and path to the SWAT+ TxtInOut folder calibration_params = dict_of_params['calibration_params'] path_to_txtinout = dict_of_params['path_to_txtinout'] # Initialize the TxtinoutReader and copy the SWAT+ project to a temporary directory reader = TxtinoutReader(path_to_txtinout) tmp_path = reader.copy_swat(dir=None) # Copy to a temporary directory reader = TxtinoutReader(tmp_path) # Run SWAT+ with the provided calibration parameters txt_in_out_result = reader.run_swat(calibration_params, show_output=False) # Initialize a new TxtinoutReader to read the results result_reader = TxtinoutReader(txt_in_out_result) """ The following steps should include: 1. Reading the simulation results. 2. Gathering observed data. 3. Calculating the error metric based on the observed and simulated data. """ # For demonstration, we return a random error metric rng = np.random.default_rng() return (rng.random(), {'test_calibration': result_reader.root_folder})

---

3. Set Up the SWAT+ Calibration Problem¶

We define the calibration problem by specifying:

• The parameters to calibrate (e.g., bm_e and harv_idx in the plants.plt file).
• The objective function (function_to_minimize).
• Additional arguments such as the path to the SWAT+ TxtInOut folder.

In [18]:

# Path to the SWAT+ TxtInOut folder
txtinout_folder = '/mnt/c/Users/joans/OneDrive/Escriptori/icra/muga_windows'

# Define the SWATProblem instance
swat_problem = SWATProblem(
    params={
        'plants.plt': ('name', [('bana', 'bm_e', 40, 50), ('bana', 'harv_idx', 0.4, 0.5)])
    },
    function_to_evaluate=function_to_minimize,
    param_arg_name='calibration_params',
    n_workers=4,
    parallelization='threads',
    debug=False,
    path_to_txtinout=txtinout_folder
)

# Path to the SWAT+ TxtInOut folder txtinout_folder = '/mnt/c/Users/joans/OneDrive/Escriptori/icra/muga_windows' # Define the SWATProblem instance swat_problem = SWATProblem( params={ 'plants.plt': ('name', [('bana', 'bm_e', 40, 50), ('bana', 'harv_idx', 0.4, 0.5)]) }, function_to_evaluate=function_to_minimize, param_arg_name='calibration_params', n_workers=4, parallelization='threads', debug=False, path_to_txtinout=txtinout_folder )

---

4. Configure the Optimization Algorithm¶

We use the CMA-ES (Covariance Matrix Adaptation Evolution Strategy) algorithm from pymoo for optimization. The algorithm is configured with:

• Initial parameter values (x0).
• Termination criteria (e.g., maximum number of evaluations).

In [19]:

# Define initial parameter values and bounds
x0 = denormalize(np.random.random(2), np.array([40, 0.4]), np.array([50, 0.5]))

# Set the number of simulations (evaluations)
n_simulations = 2

# Configure the CMA-ES algorithm
algorithm = CMAES(x0=x0, maxfevals=n_simulations)
termination = get_termination("n_eval", n_simulations)

# Define initial parameter values and bounds x0 = denormalize(np.random.random(2), np.array([40, 0.4]), np.array([50, 0.5])) # Set the number of simulations (evaluations) n_simulations = 2 # Configure the CMA-ES algorithm algorithm = CMAES(x0=x0, maxfevals=n_simulations) termination = get_termination("n_eval", n_simulations)

---

5. Run the Optimization¶

We run the optimization using the minimize_pymoo function. The results include:

• The best set of parameters (x).
• The path to the best simulation results (path).
• The error metric (error).

In [20]:

# Run the optimization
x, path, error = minimize_pymoo(
    swat_problem,
    algorithm,
    termination,
    verbose=False,
)

# Run the optimization x, path, error = minimize_pymoo( swat_problem, algorithm, termination, verbose=False, )

---

6. Analyze the Results¶

Finally, we analyze the results of the optimization:

• The best combination of parameters.
• The path to the simulation results.
• The error metric.

In [21]:

# Best combination of parameters
print("Best parameters:", x)

# Path to the best simulation results
print("Simulation results path:", path)

# Error metric
print("Error:", error)

# Best combination of parameters print("Best parameters:", x) # Path to the best simulation results print("Simulation results path:", path) # Error metric print("Error:", error)

Best parameters: [42.2548061   0.41222274]
Simulation results path: {'test_calibration': PosixPath('/tmp/tmp6gnrs7xq')}
Error: 0.2528642743696815