from shilps.powersystems import TimeConfig
from abc import ABC, abstractmethod
from typing import List, Dict, Any
from .components import PlantTemplate, DataTimeSeries
from dataclasses import dataclass
import numpy as np
import random
from deap import base, creator, tools, algorithms
[docs]
class ModelPlantDesign(ABC):
def __init__(self, plant_template) -> None:
self.plant_template = plant_template
self.plant_design = None
[docs]
@abstractmethod
def optimize(self):
pass
[docs]
class HydrogenPlantDesign(ModelPlantDesign):
[docs]
@dataclass
class SolverParams:
feas_tol: float = 1e-6
max_iter: int = 100
def __init__(self, plant_template: PlantTemplate, time_config:TimeConfig,
tsdata:DataTimeSeries) -> None:
super().__init__(plant_template)
self._nvars = 0
self._ncons = 0
self._constraint_functions = []
self._connames = np.ndarray(0, dtype=object)
self._int2idx = {}
self._idx2int = {}
self._con_int2idx = {}
self._con_idx2int = {}
self._ready_to_optimize = False
self.y_int = []
# Solver params
self.solver_params = HydrogenPlantDesign.SolverParams()
self.tsdata = tsdata
for pv_system in self.plant_template.pv_systems.values():
pv_system.solar_irradiance.tsdata = self.tsdata
# Sets
self.time_config = time_config
self.time_range = time_config.range
scenario = 1 # hardcoded scenario for now
self.periods = range(len(self.time_range))
self.set_electrolizers = list(plant_template.electrolizers.keys())
self.set_batteries = list(plant_template.batteries.keys())
self.set_pv_systems = list(plant_template.pv_systems.keys())
self.lb = np.ndarray(0)
self.ub = np.ndarray(0)
# Variables
self.add_variables("cap_el", sets=(self.set_electrolizers,),
lb=0.0)
#self.add_variables("cap_bat", sets=(self.set_batteries,),
# lb=0.0)
self.add_variables("cap_pv", sets=(self.set_pv_systems,),
lb=0.0)
#self.add_variables("pb_p", sets=(self.periods, self.set_batteries),
# lb=0.0)
#self.add_variables("pb_m", sets=(self.periods, self.set_batteries),
# lb=0.0)
#self.add_variables("sb", sets=(self.periods, self.set_batteries),
# lb=0.0)
self.add_variables("pg", sets=(self.periods, self.set_pv_systems),
lb=0.0)
self.add_variables("y", sets=(self.periods, self.set_electrolizers),
lb=0.0)
# Constraints
# Pmax value
for i in self.set_pv_systems:
for t, tt in zip(self.periods, self.time_range):
pg_idx = self._idx2int[("pg", t, i)]
f = lambda x, pg_idx=pg_idx, i=i, scenario=scenario, tt=tt : x[pg_idx] - self.plant_template.pv_systems[i].solar_irradiance[scenario,tt]
self.add_constraint(("pv_max", t, i), f)
# Gain (hardcoded)
FIRST_ELEC = self.set_electrolizers[0]
for t in self.periods:
y_int = self._idx2int[("y", t, FIRST_ELEC)]
pg_int_t = [self._idx2int[("pg", t, i)] for i in self.set_pv_systems]
f1 = lambda x, pg_int_t=pg_int_t, y_int=y_int,FIRST_ELEC=FIRST_ELEC: np.sum(self.plant_template.electrolizers[FIRST_ELEC].gain(x[pg_int_t])) - x[y_int]
self.add_constraint(("pv_gain1", t, i), f1)
f2 = lambda x: -f1(x)
self.add_constraint(("pv_gain2", t, i), f2)
self.update()
[docs]
def update(self):
self.x = np.zeros(self._nvars)
self._y_int = [self._idx2int[("y", t, i)] for t in self.periods for i in self.set_electrolizers]
[docs]
def add_constraint(self, nameindex, func):
self._constraint_functions.append(func)
self._con_idx2int[nameindex] = self._ncons
self._ncons += 1
[docs]
def get_constraint(self, nameindex):
return self._constraint_functions[self._con_idx2int[nameindex]]
[docs]
def get_constraints_names(self):
return list(self._con_idx2int.keys())
[docs]
def add_variables(self, name: str, sets: tuple, lb: float = None,
ub: float = None):
product_space = np.array(np.meshgrid(*sets)).T.reshape(-1, len(sets))
product_space = tuple(map(lambda x: tuple(int(i) for i in x), product_space))
name_indices = [(name,) + idx for idx in product_space]
self.lb = np.concatenate([self.lb, np.full(len(name_indices), lb)])
self.ub = np.concatenate([self.ub, np.full(len(name_indices), ub)])
for name_index in name_indices:
self._int2idx[self._nvars] = name_index
self._idx2int[name_index] = self._nvars
self._nvars += 1
[docs]
def compute_rhs(self, x=None):
if x is None:
x = self.x
return np.stack([f(x) for f in self._constraint_functions])
[docs]
def check_constraints(self, x):
rhs = np.stack([f(x) for f in self._constraint_functions])
return np.all(rhs <= self.solver_params.feas_tol)
[docs]
def display_variables(self):
print("{:<5s} {:<15s} {:<11s} {:<10s} {}".format(
"|Index", "|Variable", "|Lower Bound", "|Upper Bound", "|Current Value"))
for nidx, i in self._idx2int.items():
lb = self.lb[i]
lb = "{:<12s}".format("-") if lb is None else "{:<12.2f}".format(lb)
ub = self.ub[i]
ub = "{:<12s}".format("-") if ub is None else "{:<12.2f}".format(ub)
if self.x is None:
value = "-"
else:
value = self.x[i]
varname = nidx[0] + "["+ ",".join(str(x) for x in nidx[1:]) + "]"
print("{:<6s} {:<16s} {} {} {}".format(
str(i), varname, lb, ub, value))
def _evaluate(self, x):
if not self.check_constraints(x):
return - np.sum(x.clip(0.)),
return np.sum(x[self._y_int]),
[docs]
def optimize(self):
# Create the fitness and individual classes
creator.create("Fitness", base.Fitness, weights=(1.0,))
creator.create("Individual", np.ndarray, fitness=creator.Fitness)
# Initialize the toolbox
toolbox = base.Toolbox()
# Attribute generator: Generate a binary array for each individual
#toolbox.register("attr_bool", np.random.randint, 2, size=num_items)
# Attribute generator: Generate a continuous array for each individual
toolbox.register("attr_float", random.random)
# Structure initializers
toolbox.register("individual", tools.initRepeat, creator.Individual, toolbox.attr_float, n=self._nvars)
toolbox.register("population", tools.initRepeat, list, toolbox.individual)
# Evaluation function
def evaluate(x):
if not self.check_constraints(x):
return - np.sum(self.compute_rhs(x).clip(0.)),
return np.sum(x[self._y_int]),
# Register the genetic operators
toolbox.register("mate", tools.cxTwoPoint)
toolbox.register("mutate", tools.mutFlipBit, indpb=0.1)
toolbox.register("select", tools.selTournament, tournsize=3)
toolbox.register("evaluate", evaluate)
random.seed(42)
# Create an initial population
population = toolbox.population(n=300)
# Define the statistics to be gathered
stats = tools.Statistics(lambda ind: ind.fitness.values)
stats.register("avg", lambda pop: sum(f[0] for f in pop) / len(pop))
stats.register("min", min)
stats.register("max", max)
# Run the algorithm
population, log = algorithms.eaSimple(population, toolbox, cxpb=0.5, mutpb=0.2, ngen=100, stats=stats, verbose=True)
# Extract the best solution
self.x = best_individual = tools.selBest(population, 1)[0]
print("Best individual: ", best_individual)
print("Best fitness: ", best_individual.fitness.values[0])
[docs]
def validate(self):
pass
[docs]
def eval_fobj(self):
pass