Tutorial

This example is code (slightly modified) used for the Medwater project (see the technical paper here).

Python Code

Create the main script that configures and duplicates the optimization environment, then runs the optimization algorithm.

import mooflow
from eval_modflow import eval_modflow
configuration = {
     "num_parameter": 10,                    # number of optimization parameters
     "eval_function": eval_modflow,  # number of optimization parameters
     "num_threads": 12,                              # parallelisation using 12 threads
     "fitness_weights": [-1, -1],    # -1 to minimize a fitness function
     "fitnessfile_with_path": "d:/modflow_model/eval_modflow.py"
     }

opt_mocmaes = mooflow.create_mocmaes(configuration)  # initialize the optimization algorithm
mocmaes.setup_modflow_parallel(configuration)            # duplicate the model for parallelization
mocmaes.run()                                            # run the optimization

The eval_modflow.py as mentioned in fitnessfile_with_path looks like this:

# -*- coding: utf-8 -*-
from datetime import timedelta
from datetime import datetime
import mooflow as omo
import numpy

def external_modflow(parameters, threadnumber):

    #############################
    # model setup
    #############################

    # working directory for external model
    mf = "d:\\medwater_2050_rri\\Model_2020__" + str(threadnumber)
    # some supplementary files are here
    mf2 = 'd:\\medwater\\MODEL_FILES'
    # directory for output files (log output files with parameter and fitness function values go here)
    mf3 = 'd:\\medwater_2050_rri\\'
    # name of modflow executable
    exe = "d:\\medwater_2050_rri\\Model_2020__" + str(threadnumber) + "\\MODFLOW-NWT_64.exe"
    # well file: specifications for each well (x, y, z, name, wellgroup, active node, pattern class)
    wellgroupfile = 'd:\\medwater\\MODEL_FILES\\wells_list.txt'
    # extraction pattern file: 16 generalized extraction patterns
    patternfile = 'd:\\medwater\\MODEL_FILES\\pattern_16_inf.txt'
    # boundary file: the upper and lower boundaries for the patterns of each well
    boundfile = 'd:\\medwater\\MODEL_FILES\\boundaries_neu.csv'
    # pywr model file
    pywr_file = 'd:\\medwater\\MODEL_FILES\\nc_2020_des_m50__2050_rri.json'

    # create the model setup
    setup_mf = omo.classes.Omo_setup('month')
    setup_mf.modelfolder = mf
    setup_mf.usePywr = True
    setup_mf.modePywr = "total"
    setup_mf.modelname = 'WMA'
    setup_mf.starting_date = datetime(1987, 9, 15)

    # number of periods to skip for model set in
    setup_mf.setting_time = 2
    # read the temporal discretization from the dis file
    setup_mf.read_disfile()
    # sub-directory for input data
    setup_mf.inputfolder = 'ref'
    setup_mf.executable = exe


    #############################
    # Parameter class
    #############################
    Parameter = omo.classes.Parameters(parameters)

    #############################
    # registry for well groups and wells
    # several wells are grouped and share the same parameter(s)
    # which scale the extraction/infiltration pattern
    # pattern and upper/lower bounds for scaling
    # are read from a binary file
    #############################
    reg_wellgroup, reg_well = omo.read.wells_and_wellgroups(setup_mf,
       wellgroupfile, patternfile, delimiter_wg="\t", delimiter_pf="\t",
       factor=1)
    reg_wellgroup = omo.read.boundaries(boundfile, reg_wellgroup, factor=1)

   #############################
   # setup the Pywr water allocation model
   # the models configuration is stored in a json file
   # active nodes are the sectoral demands
   #############################
    Pywr_model = omo.classes.Pywr_model(setup_mf, timeseries=False)
    Pywr_model.set_pywrfile(pywr_file, mf2)
    Pywr_model.active_nodes = ['a214_IS', 'a210_IS',
            'a210_WB', 'a211_IS',
                'a211_WB', 'a212_IS',
                'a212_WB', 'a220_IS',
                'a220_WB', 'a221_WB',
                'Inf_210IS', 'Inf_211IS']

    #############################
    # ground water level sensitive locations
    # locations are areas which are set by the row/col grid numbers
    # or each area an median is calculated over all grid points-
    # finally an indicator can be selected
    # which can later be used in a fitness function
    #############################
    reg_argw = omo.classes.Register('Area_gw')
    well_file0 = omo.read.read_well_file(setup_mf, 1)
    well_locations = [
        [1, 41, 93],
        [1, 124, 71],
        [1, 310, 51],
        [1, 234, 54],
        [1, 79, 73],
        [1, 92, 106],
        [1, 172, 80],
        [2, 158, 101],
        [3, 197, 115],
        [3, 267, 103]
    ]
   well_thresholds = [9, 12, 13, 13, 13, 20, 18, 35, 295, 345]

    for i, j in enumerate(well_locations):
             for fitness_gwl in range(2, 4):
                     loc = omo.classes.Area_gw('loc' + str(i) + "_" + str(fitness_gwl),
                                 [[j[ 1], j[2]],
                                 [j[1]-1, j[2]],
                                 [j[1]+1, j[2]],
                                 [j[1], j[2]-1],
                                 [j[1], j[2]+1]])

             if fitness_gwl == 0:
                     # calculate 0.5 quantile of area
                     loc.set_area_function1(numpy.percentile, 50, axis=1)
             elif fitness_gwl == 1:
                     # max. len of all periods with ground water level below threshold
                     loc.set_area_function1(omo.funcs.max_len_ones, indicator="max_len", parameter=1)  ### BUG!
             elif fitness_gwl == 2:
                     # total number of days whith a ground water level  below the threshold
                     loc.set_area_function1(omo.funcs.days_violated, loc, setup_mf)
             else:
                     # linear trend in ground water level (times -1) if the slope is negative, else 0
                     loc.set_area_function1(omo.funcs.linear_slope, plot_flag=0)

             # the treshold in gwl can be a pattern or is constant
             loc.set_pattern([365], [well_thresholds[i]])

             # for several z layers at a location set weights for each layer
             fg = 'redlines' if fitness_gwl == 2 else 'slopes'
             loc.weights = 1
             loc.z = j[0]
             loc.fitness_group = fg
             reg_argw.add(loc)  # register the ground water level sensitive location


    #############################
    # MOCMAES algorithm is "scale sensitive"
    # all parameters are handled between 0 and 1.
    # we transform the parameters back from scaled range [0, ..., 1]
    # --> we rescale the extractions for real word usage in MODFLOW
    #############################
    for _, wg in enumerate(reg_wellgroup):
        Parameter.rescale_for_wellgroup(wg)

    #############################
    # write the model configuration for wells
    #############################
    omo.write.write_well_dinm(setup_mf, reg_well)


    #############################
    # we can now evaluate the parameter vector with MODFLOW.
    # call the local copy of MODFLOW in the subfolder ending with the right tread
    #############################
    setup_mf.run_model(silent=True)

    #############################
    # read results
    #############################

    # read the heads file and return a dictionary with the
    heads for all area_gw
    dict_heads = omo.read.read_b_heads_reg(reg_argw, setup_mf)

    #############################
    # analyse lst file
    #############################

    # read the balance section from the list file / get time series for drains and storage
    dict_lst = omo.analyse.list_file(setup_mf, reg_well)

    #############################
    # prepare Pywt model
    #############################

    # calculate the available water for input nodes, read shortages in pumping (MODFLOW NWT model required)
    dict_node_provides = Pywr_model.pywt_extr_nwt(setup_mf,
        reg_well, reg_wellgroup)

    #############################
    # prepare and run Pywt model
    #############################

    # prepare the pywr model and run it
    Pywr_model.update_nodes_in_model(1/20) # work with mean values for the 20 year perios
    Pywr_model.run_pywr()

    #############################
    # analyse the model
    #############################

    # read the results from the pywr model
    Pywr_model.read_results()

    #############################
    # compile the fitness functions
    #############################

    # the number of consecutive days with a drainage outflows (springs) under threshold
    ind_drain = omo.funcs.max_len_ones(dict_lst['drains array'] > 123456,
            indicator="max_len")

    # integrate the indicators for the heads --> no negative slope.
    # ~ for sustainability we want the trends in the groundwater levels over time
    # to be as close to zero or positive
    ff_eco = ind_drain \
        + omo.funcs.sum_neg_slopes(dict_heads, 'slopes', multiplier=1000) \
        + sum([dict_heads[x]["values"] for x in dict_heads if dict_heads[x]["fitness_group"] == "redline"])

    # the weighted sum of all deficits of sectoral demands
    # (1) deficits occur when the pumping rate was set too low (scaling of the extraction pattern)
    # or if (2) the realized extraction in MODFLOW was not sufficient or
    # (3) available alternative water sources are insufficient
    ff_shrt = sum([Pywr_model.results['Agr_210IS']['deficit_percent'],
            Pywr_model.results['Agr_211IS']['deficit_percent'],
                Pywr_model.results['Agr_212IS']['deficit_percent'],
                Pywr_model.results['Agr_214IS']['deficit_percent'],
                Pywr_model.results['Agr_220IS']['deficit_percent'],
                Pywr_model.results['Agr_210WB']['deficit_percent'],
                Pywr_model.results['Agr_211WB']['deficit_percent'],
                Pywr_model.results['Agr_212WB']['deficit_percent'],
                Pywr_model.results['Agr_220WB']['deficit_percent'],
                Pywr_model.results['Agr_221WB']['deficit_percent']]) \
                + sum([Pywr_model.results['Mun_210IS']['deficit_percent'],
                Pywr_model.results['Mun_211IS']['deficit_percent'],
                Pywr_model.results['Mun_212IS']['deficit_percent'],
                Pywr_model.results['Mun_214IS']['deficit_percent'],
                Pywr_model.results['Mun_220IS']['deficit_percent'],
                Pywr_model.results['Mun_210WB']['deficit_percent'],
                Pywr_model.results['Mun_211WB']['deficit_percent'],
                Pywr_model.results['Mun_212WB']['deficit_percent'],
                Pywr_model.results['Mun_220WB']['deficit_percent'],
                Pywr_model.results['Mun_221WB']['deficit_percent']])

    # alternative sources have different costs for providing the water'''
    ff_mon = sum([Pywr_model.results['Aschod']['provided_percent'],
            Pywr_model.results['Hadera']['provided_percent'],
                Pywr_model.results['Aschkelon']['provided_percent'],
                Pywr_model.results['Palmachim']['provided_percent'],
                Pywr_model.results['Sorek']['provided_percent']]) * 0.3 \
                + sum([Pywr_model.results['tww_210IS']['provided_percent'],
                Pywr_model.results['tww_211IS']['provided_percent'],
                Pywr_model.results['tww_212IS']['provided_percent'],
                Pywr_model.results['tww_214IS']['provided_percent'],
                Pywr_model.results['tww_220IS']['provided_percent'],
                Pywr_model.results['tww_210WB']['provided_percent'],
                Pywr_model.results['tww_211WB']['provided_percent'],
                Pywr_model.results['tww_212WB']['provided_percent'],
                Pywr_model.results['tww_220WB']['provided_percent'],
                Pywr_model.results['tww_221WB']['provided_percent'] ]) * 0.3 \
                + sum([Pywr_model.results['WFP_external_IS']['provided_percent'],
                Pywr_model.results['WFP_external_WB']['provided_percent']]) * 0.4

    # write the optimization log
    omo.utils.write_logfile(threadnumber,
        [ff_eco+dist, ff_shrt+dist, ff_mon+dist,
        ff_eco, ff_shrt, ff_mon], mf3, parameters)

    print("thread", threadnumber, "dist", dist, "eco", ff_eco, "shrt", ff_shrt, "mon", ff_mon)

    # dist is the distance to the feasible parameter space. the optimization algorihm is allowed
    # to violate this space. parameter values are clipped to space for MODFLOW.
    # the distance acts as a penalty, so the algorithm can learn the feasible parameter space.

    return (ff_eco+dist, ff_shrt+dist, ff_mon+dist)

The Pywr model in this example operates only on a total balance. Solving the water allocation with discrete time steps is possible. Since the Evaluation of the Pywr results is highly problem specific, therefore you need to write your own code to read and compute the indicators or fitness functions.

Supplementary files

Let’s look at the supplementary files.

Well file

Here is a section of the well file.

Index

column

layer

row

p_class

max

min

active_node

well_group

w1

2

245

32

0

-275.406

-2496.742

a220_IS

w220_IS

w2

2

223

36

0

1000

0

a220_IS

w220_IS

w3

2

223

37

0

-107.284

-972.6

a210_IS

w210_IS

The well file specifies several properties of every well with its identifier (Index). The location in the MODLOW model is given by layer, row and column. The generic extraction or infiltration pattern p_class is defined in the extraction pattern file. As the extraction patterns are scaled between 0.0 and 1.0, the true range for the specific well is recalculated with the min and max values. Finally, the membership of the well in a well_goup enables the sharing of optimization parameters (scaling factors) and active_node adds the extractions to an input node in the Pywr model.

The column names are:
  1. Index… well names

  2. layer… the layer (z) of the well

  3. row… the row (y) of the well

  4. column… the column(x) of the well

  5. p_class… the extraction pattern, e.g 0 maps to p_class_0 in the extraction pattern file

  6. max… min and max are used to rescale the extraction pattern to its native range (which is then scaled using the optimization parameter)

  7. min… see max

  8. active_node… the extractions of this well are added to this Pywr input node

  9. well_group… the respective well group (share parameters)

Extraction pattern file

Here is a section of the extraction pattern file.

Index

p_class_0

p_class_1

1

1.00E-05

1.01E-05

2

0.082936954

0.048898725

3

0.259970747

0.227219493

4

0.484915378

0.517017079

5

0.733234785

0.80380783

6

0.953547843

0.950135334

7

1

1

8

0.782034395

0.952269837

9

0.465708842

0.83875922

10

0.253697266

0.570840967

11

0.129050065

0.250065387

12

0.036914113

0.056042193

The extraction patterns are given for each season as end of season values. For monthly (30day) and weekly (7day) and daily time steps, the day of year is given at the single well level (see initialization of the Well class). This way, wells can have different season length and the user takes further care of grouping them.

Boundary file

Optimizing the extraction patterns means scaling the complete pattern up or down or use different scaling factors for single seasons. How much a pattern can be scaled up or down for a specific well_goup is defined with the boundary file.

wellgroup

boundary

1

2

3

4

5

6

7

8

9

10

11

12

w210_IS

lower

0.5

w210_IS

upper

2.0

w220_IS

lower

0.7

0.8

0.4

w220_IS

upper

1.3

1.8

1.4

The extraction pattern for w210_IS can be halved or doubled (from 0.5 to 2.0). For all wells group w220_IS only the seasons 5 to 7 are subject to optimization, all other season are not subject to change.