This is a simple-ish program to find the solution to optimize a watershed solution to reduce pollution for PGEAS.

pgeas_watershed_optimize.py

#!/usr/bin/python3

# Copyright (c) 2021 Sambhav Saggi
# Permission is hereby granted, free of charge, to any person
# obtaining a copy of this software and associated documentation
# files (the "Software"), to deal in the Software without
# restriction, including without limitation the rights to use,
# copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the
# Software is furnished to do so, subject to the following
# conditions:
# 
# The above copyright notice and this permission notice shall be
# included in all copies or substantial portions of the Software.
# 
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
# EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
# OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
# NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
# HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
# WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
# FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
# OTHER DEALINGS IN THE SOFTWARE.

import itertools

maxsites = 4
maxbmps = 4

# Sites are in this format: Site number: (Nitrate, Turbidity, Hazardous Substance)
sites = {
    1: (1,0,0),
    2: (1,0,0),
    3: (0,1,0),
    4: (1,0,0),
    5: (0,1,0),
    6: (0,1,0),
    7: (1,0,0),
    8: (0,1,0),
    9: (1,1,0),
    10: (0,0,1)
}

# Options are in this format: Option number: ([site number, ...], (Change in Nitrate, Change in Turbidity, Change in Hazardous substances))
options = {
    1: ([1,2,3,4,5,6,7,8,9,10],(0.75,0.5,1)),
    2: ([5,6,9],(0.75,0.5,1)),
    3: ([5],(1,0.6,1)),
    4: ([6,9],(0.75,0.5,1)),
    5: ([1,6,9],(1,0.75,1)),
    6: ([1,9],(0.8,1,1)),
    7: ([2],(0.25,1,1)),
    8: ([3,4],(0.9,0.3,1)),
    9: ([2,8],(0.8,0.6,1)),
    10: ([10],(1,1,0.5)),
    11: ([7],(0.25,1,1)),
    12: ([1],(0.95,1,1)),
    13: ([6,9],(0.5,0.3,1))
}

original = (5,5,1)

def run():
    # Create matrix
    matrix = []
    for siteset in itertools.combinations(range(0,len(sites)), maxsites):
        total = (0,0,0)
        for site in siteset:
            for option_key in options.keys():
                if (site + 1) in options[option_key][0]:
                    matrix.append((site+1,option_key))

    # Deduplicate
    matrix = list(dict.fromkeys(matrix))

    # Create up to four combinations of matrix (we're just gonna keep four for memory sake)
    print("Creating combinations of size " + str(maxbmps))
    matrix_combo = itertools.combinations_with_replacement(matrix,maxbmps)

    print("Finding best combination")

    # Find best combination
    best = sum(list(original))
    print("Score to beat: " + str(best))
    for combos in matrix_combo:
        total = 0
        for combo_index in range(0, len(combos)):
            combo = combos[combo_index]
            duplicates = [x for x, z in enumerate(combos) if z == combo]
            if len(duplicates) > 1:
                # Total is increased by a diminishing value
                for duplicate_index in range(0, len(duplicates)):
                    if combo_index > duplicates[duplicate_index]:
                        total+=sum(list(tuple(elem_1 * elem_2 for elem_1, elem_2 in zip(sites[combo[0]], tuple(x**(1/(combo_index+1)) for x in options[combo[1]][1])))))
            else:
                total+=sum(list(tuple(elem_1 * elem_2 for elem_1, elem_2 in zip(sites[combo[0]], options[combo[1]][1]))))
        if total < best:
            print("Found new best option with score " + str(total) + ": " + str(combos))
            best = total
        elif total == best:
            print("Found same option with score " + str(total) + ": " + str(combos))
        else:
            print(combos, end="     \r") # The spaces and the carriage return makes everything look pretty

if __name__ == "__main__":
    run()