GA

example.py

import time
import pandas
import numpy as np
from pywr.core import Model, Input, Output, Link, Storage
from pywr.parameters import Parameter, ParameterArrayIndexed
from pywr.recorders import NumpyArrayNodeRecorder

def create_model():
    # import flow timeseries for catchments
    flow = pandas.read_csv('Not a real flow series.csv')
    flow['Date'] = flow['Date'].apply(pandas.to_datetime)
    flow.set_index('Date', inplace=True)
    flow_parameter = ParameterArrayIndexed(flow['Flow'].values)

    model = Model(
        solver='glpk',
        start=flow.index[0],
        end=flow.index[-1],
    )

    catchment1 = Input(model, 'catchment1', min_flow=flow_parameter, max_flow=flow_parameter)
    catchment2 = Input(model, 'catchment2', min_flow=flow_parameter, max_flow=flow_parameter)

    reservoir1 = Storage(model, 'reservoir1', min_volume=3000, max_volume=20000, volume=16000)
    reservoir2 = Storage(model, 'reservoir2', min_volume=3000, max_volume=20000, volume=16000)

    # default control curve - this will be overwritten
    control_curve = np.array([0.0]*12, np.float32)

    # configure the transfer between the two reservoirs
    class TransferController(Parameter):
        def __init__(self, reservoir1, control_curve):
            super(TransferController, self).__init__()
            self.reservoir1 = reservoir1
            self.control_curve = control_curve
        def value(self, timestep, *args, **kwargs):
            percent_full = self.reservoir1._volume[0] / 20000.0  # TODO: this doesn't work with scenarios
            month = timestep.datetime.month
            control_curve = self.control_curve[month-1]
            control_curve = max(0, min(1, control_curve))

            if percent_full < control_curve:
                max_flow = 10.0
            else:
                max_flow = 0.0

            return max_flow
    controller = TransferController(reservoir1, control_curve)
    transfer = Link(model, 'transfer', max_flow=controller, cost=-500)

    demand1 = Output(model, 'demand1', max_flow=100.0, cost=-101)
    demand2 = Output(model, 'demand2', max_flow=81.0, cost=-100)

    river1 = Link(model, 'river1')
    river2 = Link(model, 'river2')

    # compensation flows from reservoirs
    compensation1 = Link(model, 'compensation1', max_flow=5.0, cost=-9999)
    compensation2 = Link(model, 'compensation2', max_flow=5.0, cost=-9998)

    terminator = Output(model, 'terminator', cost=1.0)

    catchment1.connect(reservoir1)
    catchment2.connect(reservoir2)
    reservoir1.connect(demand1)
    reservoir2.connect(demand2)
    reservoir2.connect(transfer)
    transfer.connect(reservoir1)
    reservoir1.connect(river1)
    reservoir2.connect(river2)
    river1.connect(terminator)
    river2.connect(terminator)
    compensation1.connect(terminator)
    compensation2.connect(terminator)

    recorders = {
        demand1: NumpyArrayNodeRecorder(model, demand1),
        demand2: NumpyArrayNodeRecorder(model, demand2),
        transfer: NumpyArrayNodeRecorder(model, transfer),
    }

    return model, control_curve, recorders

model, control_curve, rec = create_model()
toggle = False

def evaluate_model(genome):
    global model
    global control_curve
    global toggle
    if toggle:
        model.reset()
    toggle = True
    control_curve[:] = genome

    # run the model!
    model.run()

    # total deficit at each demand centre
    demand1_deficit = sum(100 - rec[model.nodes['demand1']].data)[0]
    demand2_deficit = sum(81 - rec[model.nodes['demand2']].data)[0]

    # financial cost of transfer
    cost = sum(rec[model.nodes['transfer']].data)[0] * 20.0

    # score used for algorithm
    score = cost + (demand1_deficit + demand2_deficit) * 1000.0
    return score

if __name__ == '__main__':
    # run the evaluation profile a few times for profiling
    control_curve = np.array([0.5]*12, np.float32)
    import time
    t0 = time.clock()
    for n in range(0, 10):
        evaluate_model(control_curve)
    print('{:.3f} seconds per evaluation'.format((time.clock()-t0)/(n+1)))

Not a real flow series.csv

run.py

import random

import numpy as np
import multiprocessing

from deap import algorithms
from deap import base
from deap import creator
from deap import tools

creator.create("FitnessMax", base.Fitness, weights=(1.0,))
creator.create("Individual", list, fitness=creator.FitnessMax)

toolbox = base.Toolbox()

def initControlCurve(cls):
    inst = cls([random.random() for x in range(0, 12)])
    return inst

toolbox.register("individual", initControlCurve, creator.Individual)
toolbox.register("population", tools.initRepeat, list, toolbox.individual)


from example import evaluate_model

def evaluate_individual(genome):
    score = evaluate_model(genome)
    return -score,

toolbox.register("evaluate", evaluate_individual)
toolbox.register("mate", tools.cxTwoPoint)
toolbox.register("mutate", tools.mutGaussian, mu=0, sigma=0.05, indpb=0.2)
toolbox.register("select", tools.selTournament, tournsize=3)

def main(num_processes=1):
    random.seed(640)

    if num_processes > 1:
        pool = multiprocessing.Pool(processes=num_processes)
        toolbox.register("map", pool.map)

    pop = toolbox.population(n=30)

    hof = tools.HallOfFame(1, similar=np.array_equal)

    stats = tools.Statistics(lambda ind: ind.fitness.values)
    stats.register("avg", np.mean)
    stats.register("std", np.std)
    stats.register("min", np.min)
    stats.register("max", np.max)

    try:
        algorithms.eaSimple(pop, toolbox, cxpb=0.5, mutpb=0.2, ngen=500, stats=stats, halloffame=hof)
    except KeyboardInterrupt:
        pass

    if num_processes > 1:    
        pool.close()

    return pop, stats, hof

if __name__ == "__main__":
    import argparse
    parser = argparse.ArgumentParser()
    parser.add_argument("-n", type=int, help="Number of processes", default=1)
    args = parser.parse_args()
    
    pop, stats, hof = main(num_processes=args.n)
    i = hof[0]
    print(i)
    
    import matplotlib.pyplot as plt
    plt.plot(np.arange(0.5, 12.5), i, '-k', linewidth=2)
    plt.xlim([0,12])
    plt.ylim([0,1])
    plt.grid(True)
    plt.show()

tasks.py

from celery import Celery, group

from example import evaluate_model
import numpy as np

import random
from deap import algorithms
from deap import base
from deap import creator
from deap import tools

app = Celery('tasks', broker='amqp://guest@localhost//')
# https://github.com/celery/celery/issues/897
app.conf.update(CELERY_RESULT_BACKEND='amqp')

@app.task
def evaluate(genome):
    return -evaluate_model(genome),

def parallel_map(function, genomes):
    function = function.func
    g = group([function.s(x) for x in genomes])
    r = g.apply_async()
    scores = r.get()
    return scores

def initControlCurve(cls):
    inst = cls([random.random() for x in range(0, 12)])
    return inst

creator.create("FitnessMax", base.Fitness, weights=(1.0,))
creator.create("Individual", list, fitness=creator.FitnessMax)

toolbox = base.Toolbox()
toolbox.register("individual", initControlCurve, creator.Individual)
toolbox.register("population", tools.initRepeat, list, toolbox.individual)
toolbox.register("evaluate", evaluate)
toolbox.register("mate", tools.cxTwoPoint)
toolbox.register("mutate", tools.mutGaussian, mu=0, sigma=0.05, indpb=0.2)
toolbox.register("select", tools.selTournament, tournsize=3)

def main():
    random.seed(640)

    toolbox.register("map", parallel_map)

    pop = toolbox.population(n=30)
    hof = tools.HallOfFame(1, similar=np.array_equal)

    stats = tools.Statistics(lambda ind: ind.fitness.values)
    stats.register("avg", np.mean)
    stats.register("std", np.std)
    stats.register("min", np.min)
    stats.register("max", np.max)

    algorithms.eaSimple(pop, toolbox, cxpb=0.5, mutpb=0.2, ngen=500, stats=stats, halloffame=hof)

if __name__ == '__main__':
    main()