diff --git a/gtep/contrib/__init__.py b/gtep/contrib/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/gtep/contrib/gtep_model_result.py b/gtep/contrib/gtep_model_result.py new file mode 100644 index 00000000..584af265 --- /dev/null +++ b/gtep/contrib/gtep_model_result.py @@ -0,0 +1,246 @@ +from gtep.gtep_model import ExpansionPlanningModel +from gtep.gtep_data import ExpansionPlanningData +from pyomo.core import TransformationFactory +from pyomo.contrib.appsi.solvers.highs import Highs +from pyomo.contrib.appsi.solvers.gurobi import Gurobi +from pyomo.environ import SolverFactory, value + +# Call dataset +data_path = "./gtep/data/5bus" +data_object = ExpansionPlanningData() +data_object.load_prescient(data_path) + +num_planning_year = 2 +num_rep_day = 2 +length_rep_day = 1 +num_commit_hour = 6 +num_dispat_min = 4 + +# Call and solve expansion planning model without reliability +mod_object = ExpansionPlanningModel( + stages=num_planning_year, + data=data_object.md, + num_reps=num_rep_day, + len_reps=length_rep_day, + num_commit=num_commit_hour, + num_dispatch=num_dispat_min, +) + + +def solve_expansion_model(mod_object): + mod_object.create_model() + TransformationFactory("gdp.bound_pretransformation").apply_to(mod_object.model) + TransformationFactory("gdp.bigm").apply_to(mod_object.model) + mod_object_opt = Gurobi() + mod_object.results = mod_object_opt.solve(mod_object.model) + + # Update critical buses and generators + # before running expansion planning model with reliability + + # ----------------------- Step 1: Identify critical nodes -------------- # + # NOTE: The nodes with largest demand are selected, but it could be changed + loads = { + bus: mod_object.model.loads[bus] + for bus in mod_object.model.loads + if mod_object.model.loads[bus] > 0 + } + sorted_buses = sorted(loads.keys(), key=lambda bus: loads[bus], reverse=True) + critical_buses = sorted_buses[ + :2 + ] # TODO: only select top 2 buses for now, but could be changed + noncritical_buses = list(set(mod_object.model.buses) - set(critical_buses)) + # print(critical_buses, noncritical_buses) + + # ----------------------- Step 2: Identify critical generators --------- # + # Calculate the power output of each generators + total_thermal_generation = {gen: 0 for gen in mod_object.model.thermalGenerators} + for inv_stage in mod_object.model.investmentStage.values(): + for rep_period in inv_stage.representativePeriod.values(): + for commit_period in rep_period.commitmentPeriod.values(): + for dispatch in commit_period.dispatchPeriod.values(): + for gen in mod_object.model.thermalGenerators: + thermalgen_value = value(dispatch.thermalGeneration[gen]) + total_thermal_generation[gen] += thermalgen_value + + total_renewable_generation = { + gen: 0 for gen in mod_object.model.renewableGenerators + } + for inv_stage in mod_object.model.investmentStage.values(): + for rep_period in inv_stage.representativePeriod.values(): + for commit_period in rep_period.commitmentPeriod.values(): + for dispatch in commit_period.dispatchPeriod.values(): + for gen in mod_object.model.renewableGenerators: + renewablegen_value = value(dispatch.renewableGeneration[gen]) + total_renewable_generation[gen] += renewablegen_value + + total_generation = {} + for gen in mod_object.model.generators: + if gen in mod_object.model.thermalGenerators: + total_generation[gen] = total_thermal_generation[gen] + elif gen in mod_object.model.renewableGenerators: + total_generation[gen] = total_renewable_generation[gen] + + # Calculate actual capacity factors of critical generators in critical nodes + average_capacity_factor = {} + for gen in mod_object.model.generators: + if gen in mod_object.model.thermalGenerators: + if total_thermal_generation[gen] != 0: + average_capacity_factor[gen] = abs( + total_thermal_generation[gen] + / ( + mod_object.model.thermalCapacity[gen] + * num_planning_year + * num_rep_day + * num_commit_hour + * num_dispat_min + ) + ) + else: + average_capacity_factor[gen] = 0 + + elif gen in mod_object.model.renewableGenerators: + if total_renewable_generation[gen] != 0: + average_capacity_factor[gen] = abs( + total_renewable_generation[gen] + / ( + max( + 0.001, + mod_object.model.renewableCapacity[gen] + * num_planning_year + * num_rep_day + * num_commit_hour + * num_dispat_min, + ) + ) + ) + else: + average_capacity_factor[gen] = 0 + # print(average_capacity_factor) + + # Select top three generators from all generators as critical generators from critical buses + # Select non-critical generators from critical buses + # TODO: the number of critical generators selected should be flexible + filtered_generation = {} + for gen, output in total_generation.items(): + bus_id = gen.split("_")[0] + bus_name = f"bus{bus_id}" + + if bus_name in critical_buses: + filtered_generation[gen] = output + + sorted_generators = sorted( + filtered_generation.items(), key=lambda x: x[1], reverse=True + ) + critical_generators = [ + gen[0] for gen in sorted_generators[:3] + ] # TODO: only select top 2 generators + non_critical_generators = [ + gen for gen in filtered_generation.keys() if gen not in critical_generators + ] + # print(critical_generators, non_critical_generators) + + # Update critical generators sets + critical_bus_gen = {bus: [] for bus in critical_buses} + critical_bus_thermalgen = {bus: [] for bus in critical_buses} + critical_bus_renewablegen = {bus: [] for bus in critical_buses} + + for gen in critical_generators: + bus_id = gen.split("_")[0] + bus_name = f"bus{bus_id}" + + if bus_name in critical_buses: + if gen in mod_object.model.generators: + critical_bus_gen[bus_name].append(gen) + + if gen in mod_object.model.thermalGenerators: + critical_bus_thermalgen[bus_name].append(gen) + + elif gen in mod_object.model.renewableGenerators: + critical_bus_renewablegen[bus_name].append(gen) + # print(critical_bus_gen, critical_bus_thermalgen, critical_bus_renewablegen) + + # Update noncritical generators sets + noncritical_bus_gen = {bus: [] for bus in critical_buses} + noncritical_bus_thermalgen = {bus: [] for bus in critical_buses} + noncritical_bus_renewablegen = {bus: [] for bus in critical_buses} + + for gen in non_critical_generators: + bus_id = gen.split("_")[0] + bus_name = f"bus{bus_id}" + + if bus_name in critical_buses: + if gen in mod_object.model.generators: + noncritical_bus_gen[bus_name].append(gen) + + if gen in mod_object.model.thermalGenerators: + noncritical_bus_thermalgen[bus_name].append(gen) + + elif gen in mod_object.model.renewableGenerators: + noncritical_bus_renewablegen[bus_name].append(gen) + # print(noncritical_bus_gen, noncritical_bus_thermalgen, noncritical_bus_renewablegen) + + # ----------------------- Step 3: Capacity failure states -------------- # + # Assign which critical generators are active in each failure state + import more_itertools + + failureStates = list( + range(1, 2 ** max(len(critical_bus_gen[bus]) for bus in critical_buses) + 1) + ) + active_generator_state = { + (bus, state): [] for bus in critical_buses for state in failureStates + } + for bus in critical_buses: + crit_bus_gen_pset = list(more_itertools.powerset(critical_bus_gen[bus]))[1:] + # print(list(crit_bus_gen_pset)) + + state_idx = 1 + for gen_set in crit_bus_gen_pset: + active_generator_state[bus, state_idx] = list(gen_set) + state_idx += 1 + + active_critical_gen = { + (bus, state): [] for bus in critical_buses for state in failureStates + } + active_critical_thermalgen = { + (bus, state): [] for bus in critical_buses for state in failureStates + } + active_critical_renewablegen = { + (bus, state): [] for bus in critical_buses for state in failureStates + } + + for bus in critical_buses: + for state in failureStates: + if (bus, state) in active_generator_state: + active_critical_gen[bus, state] = active_generator_state[bus, state] + + for (bus, state), generators in active_critical_gen.items(): + thermal_gen = [] + renewable_gen = [] + + for gen in generators: + if gen in mod_object.model.thermalGenerators: + thermal_gen.append(gen) + elif gen in mod_object.model.renewableGenerators: + renewable_gen.append(gen) + + active_critical_thermalgen[bus, state] = thermal_gen + active_critical_renewablegen[bus, state] = renewable_gen + # print(active_critical_gen, active_critical_thermalgen, active_critical_renewablegen) + + results_sets = { + "criticalBuses": critical_buses, + "noncriticalBuses": noncritical_buses, + "criticalGenerators": critical_bus_gen, + "criticalthermalGenerators": critical_bus_thermalgen, + "criticalrenewableGenerators": critical_bus_renewablegen, + "noncriticalthermalGenerators": noncritical_bus_thermalgen, + "noncriticalrenewableGenerators": noncritical_bus_renewablegen, + "activeCriticalthermalGenerators": active_critical_thermalgen, + "activeCriticalrenewableGenerators": active_critical_renewablegen, + "averageCapacityFactor": average_capacity_factor, + } + + return results_sets + + +# expansion_planning_results = solve_expansion_model(mod_object) diff --git a/gtep/contrib/gtep_reliability_algorithm.py b/gtep/contrib/gtep_reliability_algorithm.py new file mode 100644 index 00000000..21fcb0cd --- /dev/null +++ b/gtep/contrib/gtep_reliability_algorithm.py @@ -0,0 +1,149 @@ +from gtep.gtep_model import ExpansionPlanningModel +from gtep.gtep_solution import ExpansionPlanningSolution +from gtep.gtep_data import ExpansionPlanningData +from gtep.contrib.gtep_reliability_model import ExpansionPlanningModelwithReliability +from gtep.contrib.gtep_model_result import solve_expansion_model +from pyomo.core import TransformationFactory +from pyomo.contrib.appsi.solvers.highs import Highs +from pyomo.contrib.appsi.solvers.gurobi import Gurobi +from pyomo.environ import Var, Expression, SolverFactory +import csv +import more_itertools + + +# Call dataset +data_path = "./gtep/data/SanDiego" +data_object = ExpansionPlanningData() +data_object.load_prescient(data_path) + +num_planning_year = 2 +num_rep_day = 2 +length_rep_day = 1 +num_commit_hour = 6 +num_dispat_min = 4 + +# Call the expansion planning model without reliability +mod_object = ExpansionPlanningModel( + stages=num_planning_year, + data=data_object.md, + num_reps=num_rep_day, + len_reps=length_rep_day, + num_commit=num_commit_hour, + num_dispatch=num_dispat_min, +) + + +# Solve expansion planning model without reliability +# Export the results for the reliability-constrained model +expansion_model_results = solve_expansion_model(mod_object) + + +sol_object = ExpansionPlanningSolution() +sol_object.load_from_model(mod_object) +sol_object.dump_json("./gtep_pre_reliability_solution.json") + + +# Export the results of expansion planning without reliability +results_ref = [] +for var in mod_object.model.component_objects(Var): + var_name = var.name + for index in var: + values = var[index].value + results_ref.append((f"{var_name}[{index}]", values)) + +for expr in mod_object.model.component_objects(Expression): + expr_name = expr.name + for index in expr: + try: + values = expr[index]() + results_ref.append((f"{expr_name}[{index}]", values)) + except ValueError: + results_ref.append((f"{expr_name}[{index}]", None)) + +with open("optimal_variable_values_without_reliability.csv", "w", newline="") as file: + writer = csv.writer(file) + writer.writerow(["Name", "Value"]) # Header row + for row in results_ref: + writer.writerow(row) + +# Call expansion planning model with reliability +mod_object_rel = ExpansionPlanningModelwithReliability( + stages=num_planning_year, + data=data_object.md, + num_reps=num_rep_day, + len_reps=length_rep_day, + num_commit=num_commit_hour, + num_dispatch=num_dispat_min, + rel_data=expansion_model_results, +) +mod_object_rel.create_model() + + +# Calculate the probability of capacity failure state based on the probability of failure +prob_state = { + (bus, state): 1 + for bus in mod_object_rel.model.criticalBuses + for state in mod_object_rel.model.states +} +failure = { + gen: mod_object_rel.model.failureRate[gen] + for gen in mod_object_rel.model.generators +} + +for bus in mod_object_rel.model.criticalBuses: + crit_generator_sets = list( + more_itertools.powerset(mod_object_rel.model.criticalGenerators[bus]) + ) + state_idx = 1 + for cg_set in crit_generator_sets: + fail_rate = 1 + for gen in mod_object_rel.model.criticalGenerators[bus]: + if gen in cg_set: + fail_rate *= failure[gen] + else: + fail_rate *= 1 - failure[gen] + prob_state[bus, state_idx] = fail_rate + state_idx += 1 + + +# Update probability of failure of capacity failure state +for bus in mod_object_rel.model.criticalBuses: + for state in mod_object_rel.model.states: + mod_object_rel.model.prob[bus, state] = prob_state[bus, state] + + +# Transform GDP to MIP +TransformationFactory("gdp.bound_pretransformation").apply_to(mod_object_rel.model) +TransformationFactory("gdp.bigm").apply_to(mod_object_rel.model) + +# Solve expansion planning model with reliability +opt = Gurobi() +mod_object_rel.results = opt.solve(mod_object_rel.model) + + +# Export results +results_rel = [] +for var in mod_object_rel.model.component_objects(Var): + var_name = var.name + for index in var: + values = var[index].value + results_rel.append((f"{var_name}[{index}]", values)) + +for expr in mod_object_rel.model.component_objects(Expression): + expr_name = expr.name + for index in expr: + try: + values = expr[index]() + results_rel.append((f"{expr_name}[{index}]", values)) + except ValueError: + results_rel.append((f"{expr_name}[{index}]", None)) + +with open("optimal_variable_values_with_reliability.csv", "w", newline="") as file: + writer = csv.writer(file) + writer.writerow(["Name", "Value"]) # Header row + for row in results_rel: + writer.writerow(row) + +sol_object = ExpansionPlanningSolution() +sol_object.load_from_model(mod_object_rel) +sol_object.dump_json("./gtep_reliability_solution.json") diff --git a/gtep/contrib/gtep_reliability_model.py b/gtep/contrib/gtep_reliability_model.py new file mode 100644 index 00000000..c42b0e5e --- /dev/null +++ b/gtep/contrib/gtep_reliability_model.py @@ -0,0 +1,2392 @@ +# Generation and Transmission Expansion Planning with Reliability Constraints +# IDAES project +# author: Kyle Skolfield and Seolhee Cho +# date: 01/04/2024 +# Model available at http://www.optimization-online.org/DB_FILE/2017/08/6162.pdf + +from pyomo.environ import * +from pyomo.environ import units as u + +# from pyomo.gdp import * + +from egret.data.model_data import ModelData +from egret.model_library.transmission.tx_utils import scale_ModelData_to_pu +from pyomo.common.timing import TicTocTimer +from pyomo.repn.linear import LinearRepnVisitor +import json +import numpy as np + +import math + + +from math import ceil +from gtep.config_options import _get_model_config + +# Define what a USD is for pyomo units purposes +# This will be set to a base year and we will do NPV calculations +# based on automatic pyomo unit transformations +u.load_definitions_from_strings(["USD = [currency]"]) + +rng = np.random.default_rng(seed=123186) + +#################################### +########## New Work Here ########### +#################################### + +## TODO: Egret features + + +# This is only used for reporting potentially bad (i.e., large magnitude) coefficients +# and thus only when that argument is passed +class VisitorConfig(object): + def __init__(self): + self.subexpr = {} + self.var_map = {} + self.var_order = {} + + def __iter__(self): + return iter((self.subexpr, self.var_map, self.var_order)) + + +class ExpansionPlanningModelwithReliability: + """A generalized generation and transmission expansion planning model with reliability.""" + + def __init__( + self, + config=None, + stages=1, + formulation=None, + data=None, + num_reps=3, + len_reps=24, + num_commit=24, + num_dispatch=4, + rel_data=None, + ): + """Initialize generation & expansion planning model object. + + :param stages: integer number of investment periods + :param formulation: Egret stuff, to be filled + :param data: full set of model data + :param num_reps: integer number of representative periods per investment period + :param len_reps: (for now integer) length of each representative period (in hours) + :param num_commit: integer number of commitment periods per representative period + :param num_dispatch: integer number of dispatch periods per commitment period + :return: Pyomo model for full GTEP + """ + + self.stages = stages + self.formulation = formulation + self.data = data + self.num_reps = num_reps + self.len_reps = len_reps + self.num_commit = num_commit + self.num_dispatch = num_dispatch + self.config = _get_model_config() + self.timer = TicTocTimer() + self.rel_data = rel_data + + def create_model(self): + """Create concrete Pyomo model object associated with the ExpansionPlanningModel with reliability""" + + self.timer.tic("Creating GTEP Model") + m = ConcreteModel() + m.reliability_model_data = self.rel_data + + ## TODO: checks for active/built/inactive/unbuilt/etc. gen + ## NOTE: scale_ModelData_to_pu doesn't account for expansion data -- does it need to? + if self.data is None: + raise + elif type(self.data) is list: + # If self.data is a list, it is a list of data for + # representative periods + m.data_list = self.data + m.md = scale_ModelData_to_pu(self.data[0]) + else: + # If self.data is an Egret model data object, representative periods will just copy it unchanged + m.data_list = None + m.md = scale_ModelData_to_pu(self.data) + m.formulation = self.formulation + + model_set_declaration( + m, self.stages, rep_per=[i for i in range(1, self.num_reps + 1)] + ) + m.representativePeriodLength = Param( + m.representativePeriods, within=PositiveReals, default=24, units=u.hr + ) + m.numCommitmentPeriods = Param( + m.representativePeriods, + within=PositiveIntegers, + default=2, + initialize=self.num_commit, + ) + m.numDispatchPeriods = Param( + m.representativePeriods, + within=PositiveIntegers, + default=2, + initialize=self.num_dispatch, + ) + m.commitmentPeriodLength = Param(within=PositiveReals, default=1, units=u.hr) + # TODO: index by dispatch period? Certainly index by commitment period + m.dispatchPeriodLength = Param(within=PositiveReals, default=15, units=u.min) + + model_data_references(m) + model_create_investment_stages(m, self.stages) + create_objective_function(m) + + self.model = m + + ## TODO: this should handle string or i/o object for outfile + def report_model(self, outfile="pretty_model_output.txt"): + """Pretty prints Pyomo model to outfile. + + :outfile: (str, optional) _description_. Defaults to "pretty_model_output.txt". + """ + with open(outfile, "w") as outf: + self.model.pprint(ostream=outf) + + def report_large_coefficients(self, outfile, magnitude_cutoff=1e5): + """Dump very large magnitude (>= 1e5) coefficients to a json file. + + :outfile: should accept filename or open file and write there; see how we do this in pyomo elsewhere + :magnitude_cutoff: magnitude above which to report coefficients + """ + var_coef_dict = {} + for e in self.model.component_data_objects(Constraint): + cfg = VisitorConfig() + repn = LinearRepnVisitor(*cfg).walk_expression(e.body) + repn_dict = repn.linear + varname_dict = {cfg.var_map[v].name: repn_dict[v] for v in repn_dict.keys()} + var_coef_dict = dict(var_coef_dict | varname_dict) + + really_bad_var_coef_dict = { + key: value + for (key, value) in var_coef_dict.items() + if abs(value) >= magnitude_cutoff + } + really_bad_var_coef_list = sorted( + really_bad_var_coef_dict.items(), key=lambda x: x[1] + ) + with open(outfile, "w") as fil: + json.dump(really_bad_var_coef_list, fil) + + +#################################### +## Model Building Functions Below ## +#################################### + + +def add_investment_variables( + b, + investment_stage, +): + """Add variables to investment stage block. + + :param b: Investment block + :param investment_stage: Investment stage index or name + :return: None + """ + + m = b.model() + b.investmentStage = investment_stage + + # Thermal generator disjuncts (operational, installed, retired, disabled, extended) + @b.Disjunct(m.thermalGenerators) + def genOperational(disj, gen): + return + + @b.Disjunct(m.thermalGenerators) + def genInstalled(disj, gen): + return + + @b.Disjunct(m.thermalGenerators) + def genRetired(disj, gen): + return + + @b.Disjunct(m.thermalGenerators) + def genDisabled(disj, gen): + return + + @b.Disjunct(m.thermalGenerators) + def genExtended(disj, gen): + return + + @b.Disjunction(m.thermalGenerators) + def investStatus(disj, gen): + return [ + disj.genOperational[gen], + disj.genInstalled[gen], + disj.genRetired[gen], + disj.genDisabled[gen], + disj.genExtended[gen], + ] + + # Line disjuncts. For now mimicking thermal generator disjuncts, though different states may need to be defined + @b.Disjunct(m.transmission) + def branchOperational(disj, branch): + return + + @b.Disjunct(m.transmission) + def branchInstalled(disj, branch): + return + + @b.Disjunct(m.transmission) + def branchRetired(disj, branch): + return + + @b.Disjunct(m.transmission) + def branchDisabled(disj, branch): + return + + @b.Disjunct(m.transmission) + def branchExtended(disj, branch): + return + + # JSC update (done?) + # @KyleSkolfield: do we differentiate between line and transformer investments? + @b.Disjunction(m.transmission) + def branchInvestStatus(disj, branch): + return [ + disj.branchOperational[branch], + disj.branchInstalled[branch], + disj.branchRetired[branch], + disj.branchDisabled[branch], + disj.branchExtended[branch], + ] + + # Renewable generator MW values (operational, installed, retired, extended) + b.renewableOperational = Var( + m.renewableGenerators, within=NonNegativeReals, initialize=0 + ) + b.renewableInstalled = Var( + m.renewableGenerators, within=NonNegativeReals, initialize=0 + ) + b.renewableRetired = Var( + m.renewableGenerators, within=NonNegativeReals, initialize=0 + ) + b.renewableExtended = Var( + m.renewableGenerators, within=NonNegativeReals, initialize=0 + ) + + # Track and accumulate costs and penalties + b.quotaDeficit = Var(within=NonNegativeReals, initialize=0, units=u.MW) + b.operatingCostInvestment = Var(within=Reals, initialize=0, units=u.USD) + b.expansionCost = Var(within=Reals, initialize=0, units=u.USD) + b.renewableCurtailmentInvestment = Var( + within=NonNegativeReals, initialize=0, units=u.USD + ) + + Upperbounds_productions = {} + for bus in m.buses: + for state in m.states: + Upperbounds_productions[bus, state] = ( + sum( + m.renewableCapacity[renewableGen] + for renewableGen in m.renewableGenerators + ) + + sum( + m.thermalCapacity[thermalGen] for thermalGen in m.thermalGenerators + ) + + sum( + m.transmissionCapacity[line] + for line in m.transmission + if m.transmission[line]["to_bus"] == bus + ) + ) + + b.ub_prod_state = Param( + m.buses, + m.states, + domain=NonNegativeReals, + initialize=Upperbounds_productions, + doc="Upper bounds of production level from states", + ) + + def prod_state_limit(b, bus, state): + return (0, b.ub_prod_state[bus, state]) + + # Variable related to reliability + b.prod_state = Var( + m.criticalBuses, + m.states, + domain=NonNegativeReals, + bounds=prod_state_limit, + doc="Estimated production level of critical buses at failure state", + ) + + +def add_investment_constraints( + b, + investment_stage, +): + """Add standard inequalities (i.e., those not involving disjunctions) to investment stage block.""" + + m = b.model() # block for investment (the outer or largest block) + + for gen in m.thermalGenerators: + if ( + m.md.data["elements"]["generator"][gen]["in_service"] == False + and investment_stage == 1 + ): + b.genDisabled[gen].indicator_var.fix(True) + # b.genDisabled[gen].binary_indicator_var.fix(1) + elif ( + m.md.data["elements"]["generator"][gen]["in_service"] == True + and investment_stage == 1 + ): + b.genOperational[gen].indicator_var.fix(True) + # b.genInstalled[gen].binary_indicator_var.fix(1) + + for branch in m.transmission: + if ( + m.md.data["elements"]["branch"][branch]["in_service"] == False + and investment_stage == 1 + ): + b.branchDisabled[branch].indicator_var.fix(True) + # b.branchDisabled[branch].binary_indicator_var.fix(1) + elif ( + m.md.data["elements"]["branch"][branch]["in_service"] == True + and investment_stage == 1 + ): + b.branchOperational[branch].indicator_var.fix(True) + # b.branchInstalled[branch].binary_indicator_var.fix(1) + + # Planning reserve requirement constraint + ## NOTE: renewableCapacityValue is a percentage of renewableCapacity + ## TODO: renewableCapacityValue ==> renewableCapacityFactor + ## NOTE: reserveMargin is a percentage of peakLoad + ## TODO: check and re-enable with additional bounding transform before bigm + ## TODO: renewableCapacityValue... should this be time iterated? is it tech based? + ## is it site based? who can say? + """@b.Constraint() + def planning_reserve_requirement(b): + return ( + sum( + m.renewableCapacity[gen] + * m.renewableCapacityValue[gen] + * (b.renewableOperational[gen] + b.renewableInstalled[gen]) + for gen in m.renewableGenerators + ) + + sum( + m.thermalCapacity[gen] + * ( + b.genOperational[gen].indicator_var.get_associated_binary() + + b.genInstalled[gen].indicator_var.get_associated_binary() + ) + for gen in m.thermalGenerators + ) + >= (1 + m.reserveMargin[investment_stage]) * m.peakLoad[investment_stage] + )""" + + # maximum investment stage installation + ## NOTE: temporarily disabled maximum investment as default option + ## TODO: These capacities shouldn't be enabled by default since they can + ## easily cause absurd results/possibly even infeasibility. Will need to add + ## user-defined handling for this. + # @b.Constraint(m.regions) + # def maximum_thermal_investment(b, region): + # return ( + # sum( + # m.thermalCapacity[gen] + # * b.genInstalled[gen].indicator_var.get_associated_binary() + # for gen in m.thermalGenerators & m.gensAtRegion[region] + # ) + # <= b.maxThermalInvestment[region] + # ) + + # @b.Constraint(m.regions) + # def maximum_renewable_investment(b, region): + # return ( + # sum( + # m.renewableCapacity[gen] + # * b.genInstalled[gen].indicator_var.get_associated_binary() + # for gen in m.renewableGenerators & m.gensAtRegion[region] + # ) + # <= b.maxRenewableInvestment[region] + # if m.renewableGenerators & m.gensAtRegion[region] + # else Constraint.Skip + # ) + + @b.Constraint(m.criticalBuses, m.states) + def available_production_state(b, bus, state): + return b.prod_state[bus, state] == sum( + m.averageCapacityFactor[gen] + * m.renewableCapacity[gen] + * (b.renewableOperational[gen] + b.renewableInstalled[gen]) + for gen in m.noncriticalrenewableGenerators[bus] + ) + sum( + m.averageCapacityFactor[gen] + * m.thermalCapacity[gen] + * ( + b.genOperational[gen].indicator_var.get_associated_binary() + + b.genInstalled[gen].indicator_var.get_associated_binary() + ) + for gen in m.noncriticalthermalGenerators[bus] + ) + sum( + m.averageCapacityFactor[gen] + * m.renewableCapacity[gen] + * (b.renewableOperational[gen] + b.renewableInstalled[gen]) + for gen in m.activeCriticalrenewableGenerators[bus, state] + ) + sum( + m.averageCapacityFactor[gen] + * m.thermalCapacity[gen] + * ( + b.genOperational[gen].indicator_var.get_associated_binary() + + b.genInstalled[gen].indicator_var.get_associated_binary() + ) + for gen in m.activeCriticalthermalGenerators[bus, state] + ) + sum( + m.transmissionCapacity[line] + * ( + b.branchOperational[line].indicator_var.get_associated_binary() + + b.branchInstalled[line].indicator_var.get_associated_binary() + ) + for line in m.transmission + if m.transmission[line]["to_bus"] == bus + ) + + ## NOTE: The following constraints can be split into rep_per and invest_stage components if desired + + ## NOTE: Constraint (13) in the reference paper + # Minimum per-stage renewable generation requirement + @b.Constraint() + def renewable_generation_requirement(b): + renewableSurplusRepresentative = 0 + ## TODO: preprocess loads for the appropriate sum here + ed = 0 + for rep_per in b.representativePeriods: + for com_per in b.representativePeriod[rep_per].commitmentPeriods: + renewableSurplusRepresentative += ( + m.weights[rep_per] + * m.commitmentPeriodLength + * b.representativePeriod[rep_per] + .commitmentPeriod[com_per] + .renewableSurplusCommitment + ) + return ( + renewableSurplusRepresentative + b.quotaDeficit + >= m.renewableQuota[investment_stage] * ed + ) + + # Operating costs for investment period + @b.Expression() + def operatingCostInvestment(b): + operatingCostRepresentative = 0 + for rep_per in b.representativePeriods: + for com_per in b.representativePeriod[rep_per].commitmentPeriods: + operatingCostRepresentative += ( + m.weights[rep_per] + * m.commitmentPeriodLength + * b.representativePeriod[rep_per] + .commitmentPeriod[com_per] + .operatingCostCommitment + ) + return m.investmentFactor[investment_stage] * operatingCostRepresentative + + # Reliability penalty for investment period + @b.Expression() + def reliabilityPenaltyInvestment(b): + reliabilityPenaltyRepresentative = 0 + for rep_per in b.representativePeriods: + for com_per in b.representativePeriod[rep_per].commitmentPeriods: + reliabilityPenaltyRepresentative += ( + m.weights[rep_per] + * m.commitmentPeriodLength + * b.representativePeriod[rep_per] + .commitmentPeriod[com_per] + .reliabilityPenaltyCommitment + ) + return m.investmentFactor[investment_stage] * reliabilityPenaltyRepresentative + + # Investment costs for investment period + ## FIXME: investment cost definition needs to be revisited AND possibly depends on + ## data format. It is _rare_ for these values to be defined at all, let alone consistently. + @b.Constraint() + def investment_cost(b): + return b.expansionCost == m.investmentFactor[investment_stage] * ( + sum( + m.generatorInvestmentCost[gen] + * m.capitalMultiplier[gen] + * b.genInstalled[gen].indicator_var.get_associated_binary() + for gen in m.thermalGenerators + ) + + sum( + m.generatorInvestmentCost[gen] + * m.capitalMultiplier[gen] + * m.renewableCapacity[gen] + * b.renewableInstalled[gen] + for gen in m.renewableGenerators + ) + + sum( + m.generatorInvestmentCost[gen] + * m.extensionMultiplier[gen] + * b.genExtended[gen].indicator_var.get_associated_binary() + for gen in m.thermalGenerators + ) + + sum( + m.generatorInvestmentCost[gen] + * m.extensionMultiplier[gen] + * m.renewableCapacity[gen] + * b.renewableExtended[gen] + for gen in m.renewableGenerators + ) + # JSC inprog (done?) - added branch investment costs here + + sum( + m.branchInvestmentCost[branch] + * m.branchCapitalMultiplier[branch] + * b.branchInstalled[branch].indicator_var.get_associated_binary() + for branch in m.transmission + ) + + sum( + m.branchInvestmentCost[branch] + * m.branchExtensionMultiplier[branch] + * b.branchExtended[branch].indicator_var.get_associated_binary() + for branch in m.transmission + ) + ) + + # Curtailment penalties for investment period + @b.Constraint() + def renewable_curtailment_cost(b): + renewableCurtailmentRep = 0 + for rep_per in b.representativePeriods: + for com_per in b.representativePeriod[rep_per].commitmentPeriods: + renewableCurtailmentRep += ( + m.weights[rep_per] + * m.commitmentPeriodLength + * b.representativePeriod[rep_per] + .commitmentPeriod[com_per] + .renewableCurtailmentCommitment + ) + return ( + b.renewableCurtailmentInvestment + == m.investmentFactor[investment_stage] * renewableCurtailmentRep + ) + + # # Curtailment penalties for investment period + # @b.Constraint() + # def renewable_curtailment_cost(b): + # capacity_sum = 0 + # for rep_per in b.representativePeriods: + # for com_per in b.representativePeriod[rep_per].commitmentPeriods: + # for dis_per in com_per.dispatchPeriods: + # capacity_sum += ( + # b.representativePeriod[rep_per] + # .commitPeriod[com_per] + # .dispatchPeriod[dis_per] + # .capacity + # ) + + +def add_dispatch_variables( + b, + dispatch_period, +): + """Add dispatch-associated variables to representative period block.""" + + m = b.model() + c_p = b.parent_block() # dispatch -- commitment + r_p = c_p.parent_block() # dispatch -- rep + i_p = r_p.parent_block() # dispatch -- planning stage + + # Define bounds on thermal generator active generation + def thermal_generation_limits(b, thermalGen): + return (0, m.thermalCapacity[thermalGen]) + + b.thermalGeneration = Var( + m.thermalGenerators, + domain=NonNegativeReals, + bounds=thermal_generation_limits, + initialize=0, + units=u.MW * u.hr, + ) + + # Define bounds on renewable generator active generation + def renewable_generation_limits(b, renewableGen): + return (0, m.renewableCapacity[renewableGen]) + + b.renewableGeneration = Var( + m.renewableGenerators, + domain=NonNegativeReals, + bounds=renewable_generation_limits, + initialize=0, + units=u.MW * u.hr, + ) + + # Define bounds on renewable generator curtailment + def curtailment_limits(b, renewableGen): + return (0, m.renewableCapacity[renewableGen]) + + b.renewableCurtailment = Var( + m.renewableGenerators, + domain=NonNegativeReals, + bounds=curtailment_limits, + initialize=0, + units=u.MW * u.hr, + ) + + # Per generator surplus + @b.Expression(m.renewableGenerators) + def renewableGenerationSurplus(b, renewableGen): + return ( + b.renewableGeneration[renewableGen] - b.renewableCurtailment[renewableGen] + ) + + # Per generator curtailment cost + @b.Expression(m.renewableGenerators) + def renewableCurtailmentCost(b, renewableGen): + return b.renewableCurtailment[renewableGen] * m.curtailmentCost + + # Per generator cost + @b.Expression(m.thermalGenerators) + def generatorCost(b, gen): + return b.thermalGeneration[gen] * m.fuelCost[gen] + + # Load shed per bus + b.loadShed = Var(m.buses, domain=NonNegativeReals, initialize=0, units=u.MW * u.hr) + + # Per bus load shed cost + @b.Expression(m.buses) + def loadShedCost(b, bus): + return b.loadShed[bus] * m.loadShedCost + + # Variables for reliability + # These variables are indexed by investment, representative, and dispatch periods. + b.loleBuses = Var( + m.buses, + domain=NonNegativeReals, + doc="LOLE of each bus", + # TODO: unit should be added, + ) + + b.eensBuses = Var( + m.buses, + domain=NonNegativeReals, + doc="EENS of each bus", + # TODO: unit should be added, + ) + + # Reliability-related penalty + @b.Expression(m.buses) + def lole_penalty(b, bus): + return b.loleBuses[bus] * m.LOLEPenalty + + @b.Expression(m.buses) + def eens_penalty(b, bus): + return b.eensBuses[bus] * m.EENSPenalty + + # Track total dispatch values and costs + b.renewableSurplusDispatch = sum(b.renewableGenerationSurplus.values()) + + b.generationCostDispatch = sum(b.generatorCost.values()) + + b.loadShedCostDispatch = sum(b.loadShedCost.values()) + + b.curtailmentCostDispatch = sum(b.renewableCurtailmentCost.values()) + + b.eensPenaltyDispatch = sum(b.eens_penalty.values()) + + b.lolePenaltyDispatch = sum(b.lole_penalty.values()) + + b.operatingCostDispatch = ( + b.generationCostDispatch + b.loadShedCostDispatch + b.curtailmentCostDispatch + ) + + b.reliabilityPenaltyDispatch = b.eensPenaltyDispatch + b.lolePenaltyDispatch + + b.renewableCurtailmentDispatch = sum( + b.renewableCurtailment[gen] for gen in m.renewableGenerators + ) + + # Define bounds on transmission line capacity - restrictions on flow over + # uninvested lines are enforced in a disjuction below + def power_flow_limits(b, branch): + return ( + -m.transmissionCapacity[branch], + m.transmissionCapacity[branch], + ) + + # NOTE: this is an abuse of units and needs to be fixed for variable temporal resolution + b.powerFlow = Var( + m.transmission, + domain=Reals, + bounds=power_flow_limits, + initialize=0, + units=u.MW * u.hr, + ) + + @b.Disjunct(m.transmission) + def branchInUse(disj, branch): + b = disj.parent_block() + + # Voltage angle + def bus_angle_bounds(disj, bus): + return (-math.pi / 6, math.pi / 6) + + # Only create bus angle variables for the buses associated with this + # branch that is in use + disj.branch_buses = [ + bb + for bb in m.buses + if ( + m.transmission[branch]["from_bus"] == bb + or m.transmission[branch]["to_bus"] == bb + ) + ] + + disj.busAngle = Var( + disj.branch_buses, domain=Reals, initialize=0, bounds=bus_angle_bounds + ) + + # Voltage angle + def delta_bus_angle_bounds(disj, bus): + return (-math.pi / 6, math.pi / 6) + + # Rule for maximum bus angle discrepancy + def delta_bus_angle_rule(disj): + fb = m.transmission[branch]["from_bus"] + tb = m.transmission[branch]["to_bus"] + return disj.busAngle[tb] - disj.busAngle[fb] + + # @KyleSkolfield - I think this var is unused and commented it out, can we delete? + disj.deltaBusAngle = Var( + domain=Reals, bounds=delta_bus_angle_bounds, rule=delta_bus_angle_rule + ) + + ## FIXME + # @disj.Constraint() + # def max_delta_bus_angle(disj): + # return abs(disj.deltaBusAngle) <= math.pi/6 + + @disj.Constraint() + def dc_power_flow(disj): + fb = m.transmission[branch]["from_bus"] + tb = m.transmission[branch]["to_bus"] + reactance = m.md.data["elements"]["branch"][branch]["reactance"] + if m.md.data["elements"]["branch"][branch]["branch_type"] == "transformer": + reactance *= m.md.data["elements"]["branch"][branch][ + "transformer_tap_ratio" + ] + shift = m.md.data["elements"]["branch"][branch][ + "transformer_phase_shift" + ] + else: + shift = 0 + return b.powerFlow[branch] == (-1 / reactance) * ( + disj.busAngle[tb] - disj.busAngle[fb] + shift + ) + + @b.Disjunct(m.transmission) + def branchNotInUse(disj, branch): + + # JSC update (done?) Fixing power flow to 0 and not creating bus angle + # variables for branches that are not in use + @disj.Constraint() + def dc_power_flow(disj): + return b.powerFlow[branch] == 0 + + return + + # JSC update - Branches are either in-use or not. This disjunction may + # provide the basis for transmission switching in the future + @b.Disjunction(m.transmission) + def branchInUseStatus(disj, branch): + return [ + disj.branchInUse[branch], + disj.branchNotInUse[branch], + ] + + # JSC update - If a branch is in use, it must be active + # Update this when switching is implemented + @b.LogicalConstraint(m.transmission) + def must_use_active_branches(b, branch): + return b.branchInUse[branch].indicator_var.implies( + lor( + i_p.branchOperational[branch].indicator_var, + i_p.branchInstalled[branch].indicator_var, + i_p.branchExtended[branch].indicator_var, + ) + ) + + # JSC update - If a branch is not in use, it must be inactive. + # Update this when switching is implemented + @b.LogicalConstraint(m.transmission) + def cannot_use_inactive_branches(b, branch): + return b.branchNotInUse[branch].indicator_var.implies( + lor( + i_p.branchDisabled[branch].indicator_var, + i_p.branchRetired[branch].indicator_var, + ) + ) + + def lole_limit(b, bus, state): + return (0, m.dispatchPeriodLength) + + b.lole = Var( + m.criticalBuses, + m.states, + domain=NonNegativeReals, + bounds=lole_limit, + doc="LOLE at each capacity failure state", + # TODO: unit should be added, + ) + + def eens_limit(b, bus, state): + return (0, m.loads[bus]) + + b.eens = Var( + m.criticalBuses, + m.states, + domain=NonNegativeReals, + bounds=eens_limit, + doc="EENS at each capacity failure state", + # TODO: unit should be added, + ) + + # Disjunctions for reliability estimation at dispatch level + @b.Disjunct(m.criticalBuses, m.states) + def reliability_check_above(disj, bus, state): + @disj.Constraint() + def capacity_larger_than_limit(disj): + return i_p.prod_state[bus, state] >= m.loads[bus] + + @disj.Constraint() + def lole_above(disj): + return b.lole[bus, state] == 0 + + @disj.Constraint() + def eens_above(disj): + return b.eens[bus, state] == 0 + + @b.Disjunct(m.criticalBuses, m.states) + def reliability_check_below(disj, bus, state): + @disj.Constraint() + def capacity_smaller_than_limit(disj): + return i_p.prod_state[bus, state] <= m.loads[bus] + + @disj.Constraint() + def lole_below(disj): + return b.lole[bus, state] == m.dispatchPeriodLength + + @disj.Constraint() + def eens_below(disj): + return b.eens[bus, state] == m.loads[bus] - i_p.prod_state[bus, state] + + @b.Disjunction(m.criticalBuses, m.states) + def reliability_logic(disj, bus, state): + return [ + disj.reliability_check_above[bus, state], + disj.reliability_check_below[bus, state], + ] + + # Define bounds on thermal generator spinning reserve supply + def spinning_reserve_limits(b, thermalGen): + return ( + 0, + m.spinningReserveFraction[thermalGen] * m.thermalCapacity[thermalGen], + ) + + b.spinningReserve = Var( + m.thermalGenerators, + domain=NonNegativeReals, + bounds=spinning_reserve_limits, + initialize=0, + units=u.MW, + ) + + # Define bounds on thermal generator quickstart reserve supply + def quickstart_reserve_limits(b, thermalGen): + return ( + 0, + m.quickstartReserveFraction[thermalGen] * m.thermalCapacity[thermalGen], + ) + + b.quickstartReserve = Var( + m.thermalGenerators, + domain=NonNegativeReals, + bounds=quickstart_reserve_limits, + initialize=0, + units=u.MW, + ) + + +def add_dispatch_constraints(b, disp_per): + """Add dispatch-associated inequalities to representative period block.""" + m = b.model() + c_p = ( + b.parent_block() + ) # commitment block (c_p) is a parent block of current dispatch block (b) + r_p = ( + c_p.parent_block() + ) # representative block (r_p) is a parent block of commitment block (c_p) + i_p = ( + r_p.parent_block() + ) # investment block (i_p) is a parent block of representative block (r_p) + + for key in m.loads.keys(): + m.loads[key] *= max(0, rng.normal(0.5, 0.2)) + + # Energy balance constraint + @b.Constraint(m.buses) + def flow_balance(b, bus): + balance = 0 + load = m.loads.get(bus) or 0 + end_points = [ + line for line in m.transmission if m.transmission[line]["from_bus"] == bus + ] + start_points = [ + line for line in m.transmission if m.transmission[line]["to_bus"] == bus + ] + gens = [ + gen + for gen in m.generators + if m.md.data["elements"]["generator"][gen]["bus"] == bus + ] + balance -= sum(b.powerFlow[i] for i in end_points) + balance += sum(b.powerFlow[i] for i in start_points) + balance += sum(b.thermalGeneration[g] for g in gens if g in m.thermalGenerators) + balance += sum( + b.renewableGeneration[g] for g in gens if g in m.renewableGenerators + ) + balance -= load + balance += b.loadShed[bus] + return balance == 0 + + # Capacity factor constraint + # NOTE: In comparison to reference work, this is *per renewable generator* + @b.Constraint(m.renewableGenerators) + def capacity_factor(b, renewableGen): + return ( + b.renewableGeneration[renewableGen] + b.renewableCurtailment[renewableGen] + == m.renewableCapacity[renewableGen] + ) + + ## TODO: (@jkskolf) add renewableExtended to this and anywhere else + @b.Constraint(m.renewableGenerators) + def operational_renewables_only(b, renewableGen): + return ( + b.renewableGeneration[renewableGen] + <= i_p.renewableInstalled[renewableGen] + + i_p.renewableOperational[renewableGen] + ) + + # RESERVE -- total operating (spinning + quickstart) + # Total operating reserve constraint + ## NOTE: min operating reserve is a percentage of load + ## FIXME: Reserve enforcement causes infeasibility issues. We should track + ## reserve shortage in some way and find a way to penalize it -- how is this + ## done in ISOs? Is it an issue to assign this as a regional + # @b.Constraint(m.regions) + # def total_operating_reserve(b, region): + # return sum( + # b.spinningReserve[gen] + b.quickstartReserve[gen] + # for gen in m.thermalGenerators & m.gensAtRegion[region] + # ) >= m.minOperatingReserve[region] * ( + # sum( + # (m.loads.get(bus) or 0) + # for bus in m.buses + # if m.md.data["elements"]["bus"][bus]["area"] == region + # ) + # ) + + # # Total spinning reserve constraint + # @b.Constraint(m.regions) + # def total_spinning_reserve(b, region): + # return sum( + # b.spinningReserve[gen] + # for gen in m.thermalGenerators & m.gensAtRegion[region] + # ) >= m.minSpinningReserve[region] * sum( + # (m.loads.get(bus) or 0) + # for bus in m.buses + # if m.md.data["elements"]["bus"][bus]["area"] == region + # ) + + # Constraints for reliability estimation + @b.Constraint(m.criticalBuses) + def lole_critical_node(b, bus): + return b.loleBuses[bus] == sum( + m.prob[bus, state] * b.lole[bus, state] for state in m.states + ) + + @b.Constraint(m.noncriticalBuses) + def lole_noncritical_node(b, bus): + return b.loleBuses[bus] == 0 + + @b.Constraint(m.criticalBuses) + def eens_critical_node(b, bus): + return b.eensBuses[bus] == sum( + m.prob[bus, state] * b.eens[bus, state] for state in m.states + ) + + @b.Constraint(m.noncriticalBuses) + def eens_noncritical_node(b, bus): + return b.eensBuses[bus] == 0 + + +def add_commitment_variables(b, commitment_period): + """Add variables and disjuncts to commitment period block.""" + m = b.model() + r_p = b.parent_block() + i_p = r_p.parent_block() + + # Define disjunction on generator status: on/startup/shutdown/off + @b.Disjunct(m.thermalGenerators) + def genOn(disj, generator): + # operating limits + ## NOTE: Reminder: thermalMin is a percentage of thermalCapacity + b = disj.parent_block() + + # Minimum operating Limits + @disj.Constraint(b.dispatchPeriods) + def operating_limit_min(d, dispatchPeriod): + return ( + m.thermalMin[generator] + <= b.dispatchPeriod[dispatchPeriod].thermalGeneration[generator] + ) + + # Maximum operating limits + @disj.Constraint(b.dispatchPeriods) + def operating_limit_max(d, dispatchPeriod): + return ( + b.dispatchPeriod[dispatchPeriod].thermalGeneration[generator] + + b.dispatchPeriod[dispatchPeriod].spinningReserve[generator] + <= m.thermalCapacity[generator] + ) + + # Ramp up limit constraints for fully on generators + @disj.Constraint(b.dispatchPeriods, m.thermalGenerators) + def ramp_up_limits(disj, dispatchPeriod, generator): + return ( + b.dispatchPeriod[dispatchPeriod].thermalGeneration[generator] + - b.dispatchPeriod[dispatchPeriod - 1].thermalGeneration[generator] + <= m.rampUpRates[generator] + * b.dispatchPeriod[dispatchPeriod].periodLength + * m.thermalCapacity[generator] + if dispatchPeriod != 1 + else Constraint.Skip + ) + + # Ramp down limit constraints for fully on generators + @disj.Constraint(b.dispatchPeriods, m.thermalGenerators) + def ramp_down_limits(disj, dispatchPeriod, generator): + return ( + b.dispatchPeriod[dispatchPeriod - 1].thermalGeneration[generator] + - b.dispatchPeriod[dispatchPeriod].thermalGeneration[generator] + <= m.rampDownRates[generator] + * b.dispatchPeriod[dispatchPeriod].periodLength + * m.thermalCapacity[generator] + if dispatchPeriod != 1 + else Constraint.Skip + ) + + # Maximum spinning reserve constraint + ##NOTE: maxSpiningReserve is a percentage of thermalCapacity + @disj.Constraint(b.dispatchPeriods, m.thermalGenerators) + def max_spinning_reserve(disj, dispatchPeriod, generator): + return ( + b.dispatchPeriod[dispatchPeriod].spinningReserve[generator] + <= m.maxSpinningReserve[generator] * m.thermalCapacity[generator] + ) + + @b.Disjunct(m.thermalGenerators) + def genStartup(disj, generator): + b = disj.parent_block() + + # operating limits + ## NOTE: Reminder: thermalMin is a percentage of thermalCapacity + @disj.Constraint(b.dispatchPeriods) + def operating_limit_min(d, dispatchPeriod): + return 0 <= b.dispatchPeriod[dispatchPeriod].thermalGeneration[generator] + + # Maximum operating limits + @disj.Constraint(b.dispatchPeriods) + def operating_limit_max(d, dispatchPeriod): + return ( + b.dispatchPeriod[dispatchPeriod].thermalGeneration[generator] + + b.dispatchPeriod[dispatchPeriod].spinningReserve[generator] + <= m.thermalMin[generator] + ) + + # Ramp up constraints for generators starting up + ## TODO: is this max necessary? I would like to remove + @disj.Constraint(b.dispatchPeriods, m.thermalGenerators) + def ramp_up_limits(disj, dispatchPeriod, generator): + return ( + b.dispatchPeriod[dispatchPeriod].thermalGeneration[generator] + - b.dispatchPeriod[dispatchPeriod - 1].thermalGeneration[generator] + <= max( + m.thermalMin[generator], + m.rampUpRates[generator] + * b.dispatchPeriod[dispatchPeriod].periodLength, + ) + * m.thermalCapacity[generator] + if dispatchPeriod != 1 + else Constraint.Skip + ) + + @b.Disjunct(m.thermalGenerators) + def genShutdown(disj, generator): + b = disj.parent_block() + + # operating limits + ## NOTE: Reminder: thermalMin is a percentage of thermalCapacity + @disj.Constraint(b.dispatchPeriods) + def operating_limit_min(d, dispatchPeriod): + return 0 <= b.dispatchPeriod[dispatchPeriod].thermalGeneration[generator] + + # Maximum operating limits + @disj.Constraint(b.dispatchPeriods) + def operating_limit_max(d, dispatchPeriod): + return ( + b.dispatchPeriod[dispatchPeriod].thermalGeneration[generator] + + b.dispatchPeriod[dispatchPeriod].spinningReserve[generator] + <= m.thermalMin[generator] + ) + + # Ramp down constraints for generators shutting down + @disj.Constraint(b.dispatchPeriods, m.thermalGenerators) + def ramp_down_limits(disj, dispatchPeriod, generator): + return ( + b.dispatchPeriod[dispatchPeriod - 1].thermalGeneration[generator] + - b.dispatchPeriod[dispatchPeriod].thermalGeneration[generator] + <= max( + m.thermalMin[generator], + m.rampDownRates[generator] + * b.dispatchPeriod[dispatchPeriod].periodLength, + ) + * m.thermalCapacity[generator] + if dispatchPeriod != 1 + else Constraint.Skip + ) + + @b.Disjunct(m.thermalGenerators) + def genOff(disj, generator): + b = disj.parent_block() + + # operating limits + ## NOTE: Reminder: thermalMin is a percentage of thermalCapacity + @disj.Constraint(b.dispatchPeriods) + def operating_limit_max(disj, dispatchPeriod): + return b.dispatchPeriod[dispatchPeriod].thermalGeneration[generator] <= 0 + + # Maximum quickstart reserve constraint + ## NOTE: maxQuickstartReserve is a percentage of thermalCapacity + @disj.Constraint(b.dispatchPeriods, m.thermalGenerators) + def max_quickstart_reserve(disj, dispatchPeriod, generator): + return ( + b.dispatchPeriod[dispatchPeriod].quickstartReserve[generator] + <= m.maxQuickstartReserve[generator] * m.thermalCapacity[generator] + ) + + @b.Disjunction(m.thermalGenerators) + def genStatus(disj, generator): + return [ + disj.genOn[generator], + disj.genStartup[generator], + disj.genShutdown[generator], + disj.genOff[generator], + ] + + # Generators cannot be committed unless they are operational or just installed + @b.LogicalConstraint(m.thermalGenerators) + def commit_active_gens_only(b, generator): + return lor( + b.genOn[generator].indicator_var, + b.genStartup[generator].indicator_var, + b.genShutdown[generator].indicator_var, + ).implies( + lor( + i_p.genOperational[generator].indicator_var, + i_p.genInstalled[generator].indicator_var, + i_p.genExtended[generator].indicator_var, + ) + ) + + +def add_commitment_constraints( + b, + comm_per, +): + """Add commitment-associated disjunctions and constraints to representative period block.""" + m = b.model() + r_p = b.parent_block() + i_p = r_p.parent_block() + + # Define total renewable surplus/deficit for commitment block + @b.Expression() + def renewableSurplusCommitment(b): + return sum( + m.dispatchPeriodLength * b.dispatchPeriod[disp_per].renewableSurplusDispatch + for disp_per in b.dispatchPeriods + ) + + # Define total operating costs for commitment block + ## TODO: Replace this constraint with expressions using bounds transform + ## NOTE: expressions are stored in gtep_cleanup branch + ## costs considered need to be re-assessed and account for missing data + @b.Expression() + def operatingCostCommitment(b): + return ( + sum( + ## FIXME: update test objective value when this changes; ready to uncomment + # (m.dispatchPeriodLength / 60) * + b.dispatchPeriod[disp_per].operatingCostDispatch + for disp_per in b.dispatchPeriods + ) + + sum( + m.fixedOperatingCost[gen] + # * m.thermalCapacity[gen] + * ( + b.genOn[gen].indicator_var.get_associated_binary() + + b.genShutdown[gen].indicator_var.get_associated_binary() + + b.genStartup[gen].indicator_var.get_associated_binary() + ) + for gen in m.thermalGenerators + ) + ## FIXME: how do we do assign fixed operating costs to renewables; flat per location or per MW + + sum( + m.fixedOperatingCost[gen] + # * m.renewableCapacity[gen] + for gen in m.renewableGenerators + ) + + sum( + m.startupCost[gen] + * b.genStartup[gen].indicator_var.get_associated_binary() + for gen in m.thermalGenerators + ) + ) + + # Define total curtailment for commitment block + @b.Expression() + def renewableCurtailmentCommitment(b): + return sum( + b.dispatchPeriod[disp_per].renewableCurtailmentDispatch + for disp_per in b.dispatchPeriods + ) + + # Reliability penalty in commitment level + @b.Expression() + def reliabilityPenaltyCommitment(b): + return sum( + b.dispatchPeriod[disp_per].reliabilityPenaltyDispatch + for disp_per in b.dispatchPeriods + ) + + +def commitment_period_rule(b, commitment_period): + """Create commitment period block. + + :param b: commitment period block + :param commitment_period: corresponding commitment period label + """ + m = b.model() + r_p = b.parent_block() + i_p = r_p.parent_block() + + b.commitmentPeriod = commitment_period + b.dispatchPeriods = RangeSet(m.numDispatchPeriods[r_p.currentPeriod]) + b.carbonTax = Param(default=0) + b.dispatchPeriod = Block(b.dispatchPeriods) + + # update properties for this time period!!!! + if m.data_list: + m.md = m.data_list[i_p.representativePeriods.index(r_p.currentPeriod)] + + # Making an exception for cases where gens were candidates + # bc their time series reduced to single values. Will probably need to fix + # this and look at where that reduction is taking place because we need more + # than a single value if the generator is built. (Probably? Maybe there's a + # different way to handle candidate renewable data because this assumes + # knowledge of the future outputs of a candidate... could be captured by scenarios?) + # Maximum output of each renewable generator + m.renewableCapacity = { + renewableGen: ( + 0 + if type(m.md.data["elements"]["generator"][renewableGen]["p_max"]) == float + else m.md.data["elements"]["generator"][renewableGen]["p_max"]["values"][ + commitment_period - 1 + ] + ) + for renewableGen in m.renewableGenerators + } + + ## TODO: Redesign load scaling and allow nature of it as argument + # Demand at each bus + temp_scale = 3 + temp_scale = 10 + + scale_loads = True + if scale_loads: + m.loads = { + m.md.data["elements"]["load"][load_n]["bus"]: ( + temp_scale + * ( + 1 + + (temp_scale + i_p.investmentStage) / (temp_scale + len(m.stages)) + ) + ) + * m.md.data["elements"]["load"][load_n]["p_load"]["values"][ + commitment_period - 1 + ] + for load_n in m.md.data["elements"]["load"] + } + # Testing + print(m.loads) + else: + m.loads = { + m.md.data["elements"]["load"][load_n]["bus"]: m.md.data["elements"]["load"][ + load_n + ]["p_load"]["values"][commitment_period - 1] + for load_n in m.md.data["elements"]["load"] + } + + ## TODO: This feels REALLY inelegant and bad. + ## TODO: Something weird happens if I say periodLength has a unit + for period in b.dispatchPeriods: + b.dispatchPeriod[period].periodLength = Param(within=PositiveReals, default=1) + add_dispatch_variables(b.dispatchPeriod[period], period) + + add_commitment_variables(b, commitment_period) + add_commitment_constraints(b, commitment_period) + + for period in b.dispatchPeriods: + add_dispatch_constraints(b.dispatchPeriod[period], period) + + +def add_representative_period_variables(b, rep_per): + m = b.model() + i_p = b.parent_block() + + b.renewableSurplusRepresentative = Var(within=Reals, initialize=0, units=u.USD) + + +def add_representative_period_constraints(b, rep_per): + m = b.model() + i_p = b.parent_block() + + @b.LogicalConstraint(b.commitmentPeriods, m.thermalGenerators) + def consistent_commitment_shutdown(b, commitmentPeriod, thermalGen): + req_shutdown_periods = ceil( + 1 / float(m.md.data["elements"]["generator"][thermalGen]["ramp_down_rate"]) + ) + return ( + atmost( + req_shutdown_periods - 1, + [ + b.commitmentPeriod[commitmentPeriod - j - 1] + .genShutdown[thermalGen] + .indicator_var + for j in range(min(req_shutdown_periods, commitmentPeriod - 1)) + ], + ).land( + b.commitmentPeriod[commitmentPeriod - 1] + .genShutdown[thermalGen] + .indicator_var + ) + # | b.commitmentPeriod[commitmentPeriod-1].genOn.indicator_var) + .implies( + b.commitmentPeriod[commitmentPeriod] + .genShutdown[thermalGen] + .indicator_var + ) + if commitmentPeriod != 1 + else LogicalConstraint.Skip + ) + + @b.LogicalConstraint(b.commitmentPeriods, m.thermalGenerators) + def consistent_commitment_off_after_shutdown(b, commitmentPeriod, thermalGen): + req_shutdown_periods = ceil( + 1 / float(m.md.data["elements"]["generator"][thermalGen]["ramp_down_rate"]) + ) + return ( + atleast( + req_shutdown_periods, + [ + b.commitmentPeriod[commitmentPeriod - j - 1] + .genShutdown[thermalGen] + .indicator_var + for j in range(min(req_shutdown_periods, commitmentPeriod - 1)) + ], + ) + .land( + b.commitmentPeriod[commitmentPeriod - 1] + .genShutdown[thermalGen] + .indicator_var + ) + .implies( + b.commitmentPeriod[commitmentPeriod].genOff[thermalGen].indicator_var + ) + if commitmentPeriod != 1 + else LogicalConstraint.Skip + ) + + @b.LogicalConstraint(b.commitmentPeriods, m.thermalGenerators) + def consistent_commitment_startup(b, commitmentPeriod, thermalGen): + req_startup_periods = ceil( + 1 / float(m.md.data["elements"]["generator"][thermalGen]["ramp_up_rate"]) + ) + return ( + atmost( + req_startup_periods - 1, + [ + b.commitmentPeriod[commitmentPeriod - j - 1] + .genStartup[thermalGen] + .indicator_var + for j in range(min(req_startup_periods, commitmentPeriod - 1)) + ], + ).land( + b.commitmentPeriod[commitmentPeriod - 1] + .genStartup[thermalGen] + .indicator_var + ) + # | b.commitmentPeriod[commitmentPeriod-1].genOn.indicator_var) + .implies( + b.commitmentPeriod[commitmentPeriod] + .genStartup[thermalGen] + .indicator_var + ) + if commitmentPeriod != 1 + else LogicalConstraint.Skip + ) + + @b.LogicalConstraint(b.commitmentPeriods, m.thermalGenerators) + def consistent_commitment_on_after_startup(b, commitmentPeriod, thermalGen): + req_startup_periods = ceil( + 1 / float(m.md.data["elements"]["generator"][thermalGen]["ramp_up_rate"]) + ) + return ( + atleast( + req_startup_periods, + [ + b.commitmentPeriod[commitmentPeriod - j - 1] + .genStartup[thermalGen] + .indicator_var + for j in range(min(req_startup_periods, commitmentPeriod - 1)) + ], + ) + .land( + b.commitmentPeriod[commitmentPeriod - 1] + .genStartup[thermalGen] + .indicator_var + ) + .implies( + b.commitmentPeriod[commitmentPeriod].genOn[thermalGen].indicator_var + ) + if commitmentPeriod != 1 + else LogicalConstraint.Skip + ) + + @b.LogicalConstraint(b.commitmentPeriods, m.thermalGenerators) + def consistent_commitment_uptime(b, commitmentPeriod, thermalGen): + return ( + atmost( + int(m.md.data["elements"]["generator"][thermalGen]["min_up_time"]) - 1, + [ + b.commitmentPeriod[commitmentPeriod - j - 1] + .genOn[thermalGen] + .indicator_var + for j in range( + min( + int( + m.md.data["elements"]["generator"][thermalGen][ + "min_up_time" + ] + ), + commitmentPeriod - 1, + ) + ) + ], + ) + .land( + b.commitmentPeriod[commitmentPeriod - 1].genOn[thermalGen].indicator_var + ) + .implies( + b.commitmentPeriod[commitmentPeriod].genOn[thermalGen].indicator_var + ) + if commitmentPeriod + != 1 # int(m.md.data["elements"]["generator"][thermalGen]["min_up_time"])+1 + else LogicalConstraint.Skip + ) + + @b.LogicalConstraint(b.commitmentPeriods, m.thermalGenerators) + def consistent_commitment_shutdown_after_uptime(b, commitmentPeriod, thermalGen): + return ( + ( + atleast( + int(m.md.data["elements"]["generator"][thermalGen]["min_up_time"]), + [ + b.commitmentPeriod[commitmentPeriod - j - 1] + .genOn[thermalGen] + .indicator_var + for j in range( + min( + int( + m.md.data["elements"]["generator"][thermalGen][ + "min_up_time" + ] + ), + commitmentPeriod - 1, + ) + ) + ], + ).land( + b.commitmentPeriod[commitmentPeriod - 1] + .genOn[thermalGen] + .indicator_var + ) + ).implies( + b.commitmentPeriod[commitmentPeriod].genOn[thermalGen].indicator_var + | b.commitmentPeriod[commitmentPeriod] + .genShutdown[thermalGen] + .indicator_var + ) + if commitmentPeriod != 1 + else LogicalConstraint.Skip + ) + + @b.LogicalConstraint(b.commitmentPeriods, m.thermalGenerators) + def consistent_commitment_downtime(b, commitmentPeriod, thermalGen): + return ( + ( + atmost( + int(m.md.data["elements"]["generator"][thermalGen]["min_down_time"]) + - 1, + [ + b.commitmentPeriod[commitmentPeriod - j - 1] + .genOff[thermalGen] + .indicator_var + for j in range( + min( + int( + m.md.data["elements"]["generator"][thermalGen][ + "min_down_time" + ] + ), + commitmentPeriod - 1, + ) + ) + ], + ).land( + b.commitmentPeriod[commitmentPeriod - 1] + .genOff[thermalGen] + .indicator_var + ) + ).implies( + b.commitmentPeriod[commitmentPeriod].genOff[thermalGen].indicator_var + ) + if commitmentPeriod + != 1 # >= int(m.md.data["elements"]["generator"][thermalGen]["min_down_time"])+1 + else LogicalConstraint.Skip + ) + + @b.LogicalConstraint(b.commitmentPeriods, m.thermalGenerators) + def consistent_commitment_start_after_downtime(b, commitmentPeriod, thermalGen): + return ( + ( + atleast( + int( + m.md.data["elements"]["generator"][thermalGen]["min_down_time"] + ), + [ + b.commitmentPeriod[commitmentPeriod - j - 1] + .genOff[thermalGen] + .indicator_var + for j in range( + min( + int( + m.md.data["elements"]["generator"][thermalGen][ + "min_down_time" + ] + ), + commitmentPeriod - 1, + ) + ) + ], + ).land( + b.commitmentPeriod[commitmentPeriod - 1] + .genOff[thermalGen] + .indicator_var + ) + ).implies( + b.commitmentPeriod[commitmentPeriod].genOff[thermalGen].indicator_var + | b.commitmentPeriod[commitmentPeriod] + .genStartup[thermalGen] + .indicator_var + ) + if commitmentPeriod != 1 + else LogicalConstraint.Skip + ) + + +def representative_period_rule( + b, + representative_period, +): + """Create representative period block. + + :b: Representative period block + :representative_period: corresponding representative period label + """ + m = b.model() + i_s = b.parent_block() + + b.currentPeriod = representative_period + + b.commitmentPeriods = RangeSet(m.numCommitmentPeriods[representative_period]) + b.commitmentPeriod = Block(b.commitmentPeriods, rule=commitment_period_rule) + + add_representative_period_variables(b, representative_period) + add_representative_period_constraints(b, representative_period) + + +def investment_stage_rule( + b, + investment_stage, +): + """Creates investment stage block. + + :b: Investment block + :investment_stage: ID for current investment stage + """ + m = b.parent_block() + + b.representativePeriods = [ + p + for p in m.representativePeriods + # if m.representativePeriodStage[p] == investment_stage + ] + add_investment_variables(b, investment_stage) + b.representativePeriod = Block( + b.representativePeriods, rule=representative_period_rule + ) + b.maxThermalInvestment = Param(m.regions, default=1000, units=u.MW) + b.maxRenewableInvestment = Param(m.regions, default=1000, units=u.MW) + + add_investment_constraints(b, investment_stage) + + +def create_objective_function(m): + """Creates objective function. Total cost is operating cost plus + expansion cost plus penalty cost (penalties include generation deficits, + renewable quota deficits, and curtailment) + :param m: Pyomo GTEP model. + """ + if len(m.stages) > 1: + m.operatingCost = sum( + m.investmentStage[stage].operatingCostInvestment for stage in m.stages + ) + m.expansionCost = sum( + m.investmentStage[stage].expansionCost for stage in m.stages + ) + m.penaltyCost = sum( + m.deficitPenalty[stage] + * m.investmentFactor[stage] + * m.investmentStage[stage].quotaDeficit + + m.investmentStage[stage].renewableCurtailmentInvestment + for stage in m.stages + ) + + m.reliabilitypenaltyCost = sum( + m.investmentStage[stage].reliabilityPenaltyInvestment for stage in m.stages + ) + + @m.Objective() + def total_cost_objective_rule(m): + if len(m.stages) > 1: + return ( + m.operatingCost + + m.expansionCost + + m.penaltyCost + + m.reliabilitypenaltyCost + ) + else: + return ( + m.investmentStage[1].operatingCostInvestment + + m.investmentStage[1].expansionCost + + m.deficitPenalty[1] + * m.investmentFactor[1] + * m.investmentStage[1].quotaDeficit + + m.investmentStage[1].renewableCurtailmentInvestment + + m.investmentStage[1].reliabilityPenaltyInvestment + ) + + +def model_data_references(m): + """Creates and labels data for GTEP model; ties input data + to model directly. + :param m: Pyomo model object + """ + + # Maximum output of each thermal generator + m.thermalCapacity = { + thermalGen: m.md.data["elements"]["generator"][thermalGen]["p_max"] + for thermalGen in m.thermalGenerators + } + + # Lifetime of each generator; needs units + m.lifetimes = { + gen: m.md.data["elements"]["generator"][gen]["lifetime"] for gen in m.generators + } + + # Minimum output of each thermal generator + m.thermalMin = { + thermalGen: m.md.data["elements"]["generator"][thermalGen]["p_min"] + for thermalGen in m.thermalGenerators + } + + # Maximum output of each renewable generator + m.renewableCapacity = { + renewableGen: ( + 0 + if type(m.md.data["elements"]["generator"][renewableGen]["p_max"]) == float + else max( + m.md.data["elements"]["generator"][renewableGen]["p_max"]["values"] + ) + ) + for renewableGen in m.renewableGenerators + } + + # A fraction of renewableCapacity representing fraction of capacity + # that can be reliably counted toward planning reserve requirement + # TODO: WHAT HAVE I DONE HERE I HATE IT and JSC made it worse... + m.renewableCapacityValue = { + renewableGen: ( + 0 + if type(m.md.data["elements"]["generator"][renewableGen]["p_max"]) == float + else min( + m.md.data["elements"]["generator"][renewableGen]["p_max"]["values"] + ) + / max(1, m.renewableCapacity[renewableGen]) + ) + for renewableGen in m.renewableGenerators + } + + # Long term thermal rating of each transmission line + m.transmissionCapacity = { + transmissionLine: m.md.data["elements"]["branch"][transmissionLine][ + "rating_long_term" + ] + for transmissionLine in m.transmission.keys() + } + + # Maximum fraction of a thermal generator's maximum output that can be + # supplied as spinning reserve + m.spinningReserveFraction = { + thermalGen: m.md.data["elements"]["generator"][thermalGen][ + "spinning_reserve_frac" + ] + for thermalGen in m.thermalGenerators + } + + # Maximum fraction of a thermal generator's maximum output that can be + # supplied as quickstart reserve + m.quickstartReserveFraction = { + thermalGen: m.md.data["elements"]["generator"][thermalGen][ + "quickstart_reserve_frac" + ] + for thermalGen in m.thermalGenerators + } + + # Demand at each bus + m.loads = { + m.md.data["elements"]["load"][load_n]["bus"]: m.md.data["elements"]["load"][ + load_n + ]["p_load"] + for load_n in m.md.data["elements"]["load"] + } + + ## NOTE: lazy fixing for dc_branch and branch... but should be an ok lazy fix + # Per-distance-unit multiplicative loss rate for each transmission line + m.lossRate = { + branch: (m.md.data["elements"]["branch"][branch].get("loss_rate") or 0) + for branch in m.transmission + } + + ## NOTE: lazy fixing for dc_branch and branch... but should be an ok lazy fix + # Distance between terminal buses for each transmission line + m.distance = { + branch: (m.md.data["elements"]["branch"][branch].get("distance") or 0) + for branch in m.transmission + } + + # JSC TODO: Add cost of investment in each new branch to input data. Currently + # selected 0 to ensure investments will be selected if needed + m.branchInvestmentCost = { + branch: (m.md.data["elements"]["branch"][branch].get("capital_cost") or 0) + for branch in m.transmission + } + + # JSC TODO: Add branch capital multiplier to input data. + m.branchCapitalMultiplier = { + branch: (m.md.data["elements"]["branch"][branch].get("capital_multiplier") or 1) + for branch in m.transmission + } + + # Cost of life extension for each generator, expressed as a fraction of initial investment cost + m.branchExtensionMultiplier = { + branch: ( + m.md.data["elements"]["branch"][branch].get("extension_multiplier") or 1 + ) + for branch in m.transmission + } + + ## TODO: These should go into each stage -- check where these values should come from + m.peakLoad = Param(m.stages, default=0, units=u.MW) + m.reserveMargin = Param(m.stages, default=0, units=u.MW) + m.renewableQuota = Param(m.stages, default=0, units=u.MW) + m.weights = Param(m.representativePeriods, default=1) + m.investmentFactor = Param(m.stages, default=1, mutable=True) + ## NOTE: Lazy approx for NPV + ## TODO: don't lazily approx NPV, add it into unit handling and calculate from actual time frames + for stage in m.stages: + m.investmentFactor[stage] *= 1 / ((1.04) ** (5 * stage)) + m.fixedOperatingCost = Param(m.generators, default=1, units=u.USD) + m.deficitPenalty = Param(m.stages, default=1, units=u.USD / u.MW) + + # Amount of fuel required to be consumed for startup process for each generator + m.startFuel = { + gen: m.md.data["elements"]["generator"][gen]["start_fuel"] + for gen in m.generators + } + + # Cost per unit of fuel at each generator + if "RTS-GMLC" in m.md.data["system"]["name"]: + m.fuelCost = { + gen: m.md.data["elements"]["generator"][gen]["fuel_cost"] + for gen in m.thermalGenerators + } + else: + m.fuelCost = { + gen: m.md.data["elements"]["generator"][gen]["p_cost"]["values"][1] + for gen in m.thermalGenerators + } + + # Cost per MW of curtailed renewable energy + # NOTE: what should this be valued at? This being both curtailment and load shed. + # TODO: update valuations + m.curtailmentCost = 2 * max(m.fuelCost.values()) + m.loadShedCost = 1000 * m.curtailmentCost + + # Full lifecycle CO_2 emission factor for each generator + m.emissionsFactor = { + gen: m.md.data["elements"]["generator"][gen]["emissions_factor"] + for gen in m.generators + } + + # Flat startup cost for each generator + if "RTS-GMLC" in m.md.data["system"]["name"]: + m.startupCost = { + gen: m.md.data["elements"]["generator"][gen]["non_fuel_startup_cost"] + for gen in m.thermalGenerators + } + else: + m.startupCost = { + gen: m.md.data["elements"]["generator"][gen]["startup_cost"] + for gen in m.generators + } + + # (Arbitrary) multiplier for new generator investments corresponds to depreciation schedules + # for individual technologies; higher values are indicative of slow depreciation + m.capitalMultiplier = { + gen: m.md.data["elements"]["generator"][gen]["capital_multiplier"] + for gen in m.generators + } + + # Cost of life extension for each generator, expressed as a fraction of initial investment cost + m.extensionMultiplier = { + gen: m.md.data["elements"]["generator"][gen]["extension_multiplier"] + for gen in m.generators + } + + # Cost of investment in each new generator + m.generatorInvestmentCost = { + # gen: m.md.data["elements"]["generator"][gen]["investment_cost"] + # for gen in m.generators] + gen: 0 + for gen in m.generators + } + + # Minimum operating reserve, expressed as a fraction of load within a region + m.minOperatingReserve = { + region: m.md.data["system"]["min_operating_reserve"] for region in m.regions + } + + # Minimum spinning reserve, expressed as a fraction of load within a region + m.minSpinningReserve = { + region: m.md.data["system"]["min_spinning_reserve"] for region in m.regions + } + + # Maximum spinning reserve available for each generator; expressed as a fraction + # maximum generator output + m.maxSpinningReserve = { + gen: m.md.data["elements"]["generator"][gen]["max_spinning_reserve"] + for gen in m.thermalGenerators + } + + # Maximum quickstart reserve available for each generator; expressed as a fraction + # maximum generator output + m.maxQuickstartReserve = { + gen: m.md.data["elements"]["generator"][gen]["max_quickstart_reserve"] + for gen in m.thermalGenerators + } + + # Ramp up rates for each generator; expressed as a fraction of maximum generator output + m.rampUpRates = { + thermalGen: m.md.data["elements"]["generator"][thermalGen]["ramp_up_rate"] + for thermalGen in m.thermalGenerators + } + + # Ramp down rates for each generator; expressed as a fraction of maximum generator output + m.rampDownRates = { + thermalGen: m.md.data["elements"]["generator"][thermalGen]["ramp_down_rate"] + for thermalGen in m.thermalGenerators + } + + # Matching for each generator to the region containing the bus at which the generator + # is located + m.gensAtRegion = { + region: gen + for region in m.regions + for gen in m.generators + if m.md.data["elements"]["bus"][m.md.data["elements"]["generator"][gen]["bus"]][ + "area" + ] + == region + } + + # Parameters for reliability + m.failureRate = { + gen: m.md.data["elements"]["generator"][gen]["failure_rate"] + for gen in m.generators + } + + m.prob = Param(m.criticalBuses, m.states, initialize=1, mutable=True) + + m.EENSPenalty = Param(default=5) # CAISO's VOLL + m.LOLEPenalty = Param(default=5) + # (cho) NOTE: LOLE is now penalized in the objective function, but it could be added as a constraint + # TODO: Correct value for LOLE penlaty should be collected + m.averageCapacityFactor = Param( + m.generators, + initialize=m.reliability_model_data["averageCapacityFactor"], + mutable=True, + ) + # (cho) TODO: This average capacity factor should be updated depending on type of generators + + # m.CriticalgeneratorProduction = Param(m.generators, default=0, mutable=True) + + +def model_set_declaration(m, stages, rep_per=["a", "b"], com_per=2, dis_per=2): + """ + Creates Pyomo Sets necessary (convenient) for solving the GTEP model. + + :m: Pyomo model object + :stages: Number of stages in investment horizon + """ + + # Demand at each bus + m.loads = { + m.md.data["elements"]["load"][load_n]["bus"]: m.md.data["elements"]["load"][ + load_n + ]["p_load"] + for load_n in m.md.data["elements"]["load"] + } + # Call load data first as it is required to define some sets + + m.buses = Set( + initialize=m.md.data["elements"]["bus"].keys(), doc="Individual buses" + ) + + m.regions = Set( + initialize=( + m.md.data["elements"]["bus"][bus]["area"] + for bus in m.md.data["elements"]["bus"] + ), + doc="Regions / clusters of buses", + ) + + ## TODO: Right now, this means that branches can only be specified entirely as standard + ## or as dc ... not mix-and-match + if len(m.md.data["elements"]["branch"]) == 0: + m.md.data["elements"]["branch"] = m.md.data["elements"]["dc_branch"] + + m.transmission = { + branch: { + "from_bus": m.md.data["elements"]["branch"][branch]["from_bus"], + "to_bus": m.md.data["elements"]["branch"][branch]["to_bus"], + "reactance": m.md.data["elements"]["branch"][branch]["reactance"], + } + for branch in m.md.data["elements"]["branch"] + } + + m.generators = Set( + initialize=m.md.data["elements"]["generator"].keys(), + doc="All generators", + ) + + m.thermalGenerators = Set( + within=m.generators, + initialize=( + gen + for gen in m.generators + if m.md.data["elements"]["generator"][gen]["generator_type"] == "thermal" + ), + doc="Thermal generators; subset of all generators", + ) + + m.renewableGenerators = Set( + within=m.generators, + initialize=( + gen + for gen in m.generators + if m.md.data["elements"]["generator"][gen]["generator_type"] == "renewable" + ), + doc="Renewable generators; subset of all generators", + ) + + ## NOTE: will want to cover baseline generator types in IDAES + + if m.md.data["elements"].get("storage"): + m.storage = Set( + initialize=(batt for batt in m.md.data["elements"]["storage"]), + doc="Potential storage units", + ) + + ## TODO: make sure time units are both definable and consistent without being forced + + m.stages = RangeSet(stages, doc="Set of planning periods") + + m.representativePeriods = Set( + initialize=rep_per, + doc="Set of representative periods for each planning period", + ) + + # Sets for reliability + # TODO: should be flexible depending on critical nodes and generators + + m.criticalBuses = Set( + within=m.buses, + initialize=m.reliability_model_data["criticalBuses"], + doc="Critical buses; subset of buses", + ) + + m.noncriticalBuses = Set( + within=m.buses, + initialize=m.reliability_model_data["noncriticalBuses"], + doc="Non-critical buses; subset of buses, remaining sets", + ) + + m.criticalGenerators = Set( + m.criticalBuses, + within=m.generators, + initialize=m.reliability_model_data["criticalGenerators"], + doc="Critical generators; subset of all generators, initially empty", + ) + + m.criticalthermalGenerators = Set( + m.criticalBuses, + within=m.thermalGenerators, + initialize=m.reliability_model_data["criticalthermalGenerators"], + doc="Critical thermal generators; subset of all generators, initially empty", + ) + + m.criticalrenewableGenerators = Set( + m.criticalBuses, + within=m.renewableGenerators, + initialize=m.reliability_model_data["criticalrenewableGenerators"], + doc="Critical renewable generators; subset of all generators, initially empty", + ) + + m.noncriticalthermalGenerators = Set( + m.criticalBuses, + within=m.thermalGenerators, + initialize=m.reliability_model_data["noncriticalthermalGenerators"], + doc="Non-critical thermal generators; subset of all generators, initially empty", + ) + + m.noncriticalrenewableGenerators = Set( + m.criticalBuses, + within=m.renewableGenerators, + initialize=m.reliability_model_data["noncriticalrenewableGenerators"], + doc="Non-critical renewable generators; subset of all generators, initially empty", + ) + + ## FIXME: lol + m.criticalGenerators.pprint() + print(list(m.criticalGenerators[bus] for bus in m.criticalBuses)) + print(list(m.criticalGenerators[bus].data() for bus in m.criticalBuses)) + failure_state = list( + range( + 1, + 2 ** (max(len(m.criticalGenerators[bus].data()) for bus in m.criticalBuses)) + + 1, + ) + ) + m.states = Set(initialize=failure_state, doc="capacity failure states") + + m.activeCriticalthermalGenerators = Set( + m.criticalBuses, + m.states, + within=m.thermalGenerators, + initialize=m.reliability_model_data["activeCriticalthermalGenerators"], + doc="Active critical thermal generators in the state, initially empty", + ) + + m.activeCriticalrenewableGenerators = Set( + m.criticalBuses, + m.states, + within=m.renewableGenerators, + initialize=m.reliability_model_data["activeCriticalrenewableGenerators"], + doc="Active critical renewable generators in the state, initially empty", + ) + + +def model_create_investment_stages(m, stages): + """Creates investment blocks and linking constraints for GTEP model. + Largely manages retirements and links operational units in a given investment stage + to operational + installed - retired in the previous investment stage. + + :m: Pyomo model object + :stages: Number of investment stages in planning horizon + """ + + m.investmentStage = Block(m.stages, rule=investment_stage_rule) + + # Retirement/extension relationships over investment periods -- C&P'd + # from the paper. These are okay. + # if len(m.stages) > 1: + + # @m.Constraint(m.stages, m.thermalGenerators) + # def gen_retirement(m, stage, gen): + # return sum( + # m.investmentStage[t_2] + # .genInstalled[gen] + # .indicator_var.get_associated_binary() + # for t_2 in m.stages + # if t_2 <= stage - m.lifetimes[gen] + # ) <= sum( + # m.investmentStage[t_1] + # .genRetired[gen] + # .indicator_var.get_associated_binary() + # + m.investmentStage[t_1] + # .genExtended[gen] + # .indicator_var.get_associated_binary() + # for t_1 in m.stages + # if t_1 <= stage + # ) + + # Linking generator investment status constraints + @m.Constraint(m.stages, m.thermalGenerators) + def gen_stats_link(m, stage, gen): + return ( + m.investmentStage[stage] + .genOperational[gen] + .indicator_var.get_associated_binary() + == m.investmentStage[stage - 1] + .genOperational[gen] + .indicator_var.get_associated_binary() + + m.investmentStage[stage - 1] + .genInstalled[gen] + .indicator_var.get_associated_binary() + - m.investmentStage[stage - 1] + .genRetired[gen] + .indicator_var.get_associated_binary() + if stage != 1 + else Constraint.Skip + ) + + # Renewable generation (in MW) retirement relationships + if len(m.stages) > 1: + + @m.Constraint(m.stages, m.renewableGenerators) + def renewable_retirement(m, stage, gen): + return sum( + m.investmentStage[t_2].renewableInstalled[gen] + for t_2 in m.stages + if t_2 <= stage - m.lifetimes[gen] + ) <= sum( + m.investmentStage[t_1].renewableRetired[gen] + + m.investmentStage[t_1].renewableExtended[gen] + for t_1 in m.stages + if t_1 <= stage + ) + + # Total renewable generation (in MW) operational at a given stage + # is equal to what was operational and/or installed in the previous stage + # less what was retired in the previous stage + @m.Constraint(m.stages, m.renewableGenerators) + def renewable_stats_link(m, stage, gen): + return ( + m.investmentStage[stage].renewableOperational[gen] + == m.investmentStage[stage - 1].renewableOperational[gen] + + m.investmentStage[stage - 1].renewableInstalled[gen] + - m.investmentStage[stage - 1].renewableRetired[gen] + if stage != 1 + else Constraint.Skip + ) + + # If a gen is online at time t, it must have been online or installed at time t-1 + @m.LogicalConstraint(m.stages, m.thermalGenerators) + def consistent_operation(m, stage, gen): + return ( + m.investmentStage[stage] + .genOperational[gen] + .indicator_var.implies( + m.investmentStage[stage - 1].genOperational[gen].indicator_var + | m.investmentStage[stage - 1].genInstalled[gen].indicator_var + ) + if stage != 1 + else LogicalConstraint.Skip + ) + + # If a gen is online at time t, it must be online, extended, or retired at time t+1 + @m.LogicalConstraint(m.stages, m.thermalGenerators) + def consistent_operation_future(m, stage, gen): + return ( + m.investmentStage[stage - 1] + .genOperational[gen] + .indicator_var.implies( + m.investmentStage[stage].genOperational[gen].indicator_var + | m.investmentStage[stage].genExtended[gen].indicator_var + | m.investmentStage[stage].genRetired[gen].indicator_var + ) + if stage != 1 + else LogicalConstraint.Skip + ) + + # Retirement in period t-1 implies disabled in period t + @m.LogicalConstraint(m.stages, m.thermalGenerators) + def full_retirement(m, stage, gen): + return ( + m.investmentStage[stage - 1] + .genRetired[gen] + .indicator_var.implies( + m.investmentStage[stage].genDisabled[gen].indicator_var + ) + if stage != 1 + else LogicalConstraint.Skip + ) + + # If a gen is disabled at time t-1, it must stay disabled at time t + @m.LogicalConstraint(m.stages, m.thermalGenerators) + def consistent_disabled(m, stage, gen): + return ( + m.investmentStage[stage - 1] + .genDisabled[gen] + .indicator_var.implies( + m.investmentStage[stage].genDisabled[gen].indicator_var + | m.investmentStage[stage].genInstalled[gen].indicator_var + ) + if stage != 1 + else LogicalConstraint.Skip + ) + + # If a gen is extended at time t-1, it must stay extended or be retired at time t + @m.LogicalConstraint(m.stages, m.thermalGenerators) + def consistent_extended(m, stage, gen): + return ( + m.investmentStage[stage - 1] + .genExtended[gen] + .indicator_var.implies( + m.investmentStage[stage].genExtended[gen].indicator_var + | m.investmentStage[stage].genRetired[gen].indicator_var + ) + if stage != 1 + else LogicalConstraint.Skip + ) + + # Installation in period t-1 implies operational in period t + @m.LogicalConstraint(m.stages, m.thermalGenerators) + def full_investment(m, stage, gen): + return ( + m.investmentStage[stage - 1] + .genInstalled[gen] + .indicator_var.implies( + m.investmentStage[stage].genOperational[gen].indicator_var + ) + if stage != 1 + else LogicalConstraint.Skip + ) + + # If a branch is online at time t, it must have been online or installed at time t-1 + @m.LogicalConstraint(m.stages, m.transmission) + def consistent_branch_operation(m, stage, branch): + return ( + m.investmentStage[stage] + .branchOperational[branch] + .indicator_var.implies( + m.investmentStage[stage - 1].branchOperational[branch].indicator_var + | m.investmentStage[stage - 1].branchInstalled[branch].indicator_var + ) + if stage != 1 + else LogicalConstraint.Skip + ) + + # If a branch is online at time t, it must be online, extended, or retired at time t+1 + @m.LogicalConstraint(m.stages, m.transmission) + def consistent_branch_operation_future(m, stage, branch): + return ( + m.investmentStage[stage - 1] + .branchOperational[branch] + .indicator_var.implies( + m.investmentStage[stage].branchOperational[branch].indicator_var + | m.investmentStage[stage].branchExtended[branch].indicator_var + | m.investmentStage[stage].branchRetired[branch].indicator_var + ) + if stage != 1 + else LogicalConstraint.Skip + ) + + # Retirement in period t-1 implies disabled in period t + @m.LogicalConstraint(m.stages, m.transmission) + def full_branch_retirement(m, stage, branch): + return ( + m.investmentStage[stage - 1] + .branchRetired[branch] + .indicator_var.implies( + m.investmentStage[stage].branchDisabled[branch].indicator_var + ) + if stage != 1 + else LogicalConstraint.Skip + ) + + # If a branch is disabled at time t-1, it must stay disabled or be installed at time t + @m.LogicalConstraint(m.stages, m.transmission) + def consistent_branch_disabled(m, stage, branch): + return ( + m.investmentStage[stage - 1] + .branchDisabled[branch] + .indicator_var.implies( + m.investmentStage[stage].branchDisabled[branch].indicator_var + | m.investmentStage[stage].branchInstalled[branch].indicator_var + ) + if stage != 1 + else LogicalConstraint.Skip + ) + + # If a branch is extended at time t-1, it must stay extended or be retired at time t + @m.LogicalConstraint(m.stages, m.transmission) + def consistent_branch_extended(m, stage, branch): + return ( + m.investmentStage[stage - 1] + .branchExtended[branch] + .indicator_var.implies( + m.investmentStage[stage].branchExtended[branch].indicator_var + | m.investmentStage[stage].branchRetired[branch].indicator_var + ) + if stage != 1 + else LogicalConstraint.Skip + ) + + # Installation in period t-1 implies operational in period t + @m.LogicalConstraint(m.stages, m.transmission) + def full_branch_investment(m, stage, branch): + return ( + m.investmentStage[stage - 1] + .branchInstalled[branch] + .indicator_var.implies( + m.investmentStage[stage].branchOperational[branch].indicator_var + ) + if stage != 1 + else LogicalConstraint.Skip + ) diff --git a/gtep/contrib/gtep_solution.py b/gtep/contrib/gtep_solution.py new file mode 100644 index 00000000..4a938a5e --- /dev/null +++ b/gtep/contrib/gtep_solution.py @@ -0,0 +1,1312 @@ +# Generation and Transmission Expansion Planning +# IDAES project +# author: Kyle Skolfield, Thom R. Edwards +# date: 01/04/2024 +# Model available at http://www.optimization-online.org/DB_FILE/2017/08/6162.pdf + +from pyomo.environ import * +from pyomo.environ import units as u +from gtep.gtep_model import ExpansionPlanningModel +import logging + +import json +from pathlib import Path + +import matplotlib.pyplot as plt +from matplotlib.ticker import MaxNLocator +import networkx as nx +import pandas as pd +import numpy as np +import re + + +from matplotlib.patches import Rectangle, RegularPolygon, PathPatch +from matplotlib.collections import PatchCollection +import matplotlib.cm as cm +from matplotlib.transforms import Affine2D +from matplotlib.colors import Normalize +import matplotlib.path as mpath + +logger = logging.getLogger(__name__) + + +# [TODO] inject units into plots +class ExpansionPlanningSolution: + def __init__(self): + # PopPop the Power Optimization Possum says ౿ᓕ ̤Ꜥ·> --- "eat trash, heck __init__, hold me" + pass + + def load_from_file(self): + pass + + def load_from_model(self, gtep_model): + if type(gtep_model) is not ExpansionPlanningModel: + logger.warning( + f"Solutions must be loaded from ExpansionPlanningModel objects, not %s" + % type(gtep_model) + ) + raise ValueError + if gtep_model.results is None: + raise ValueError( + "ExpansionPlanningSolution objects loaded from model must have a results component." + ) + self.results = gtep_model.results # Highs results object + self.stages = gtep_model.stages # int + self.formulation = gtep_model.formulation # None (???) + self.data = gtep_model.data # ModelData object + self.num_reps = gtep_model.num_reps # int + self.len_reps = gtep_model.len_reps # int + self.num_commit = gtep_model.num_commit # int + self.num_dispatch = gtep_model.num_dispatch # int + + self.expressions = { + expr.name: value(expr) + for expr in gtep_model.model.component_data_objects(Expression) + if ("Commitment" in expr.name) or ("Investment" in expr.name) + } + + self.variables = { + expr.name: value(expr) + for expr in gtep_model.model.component_data_objects(Expression) + if ("Commitment" in expr.name) or ("Investment" in expr.name) + } + + def import_data_object(self, data_obj): + self.data = data_obj.md + + def read_json(self, filepath): + # read a json file and recover a solution primals + json_filepath = Path(filepath) + with open(json_filepath, "r") as fobj: + json_read = json.loads(fobj.read()) + self.primals_tree = json_read["results"]["primals_tree"] + + def dump_json(self, filename="./gtep_solution_jscTest.json"): + # def dump_json(self, filename="./gtep_solution.json"): + + dump_filepath = Path(filename) + with open(dump_filepath, "w") as fobj: + json.dump(self._to_dict(), fobj) + + def _to_dict(self) -> dict: + + results_dict = { + "solution_loader": self.results.solution_loader, # object + "termination_condition": self.results.termination_condition, # object + "best_feasible_objective": self.results.best_feasible_objective, + "best_objective_bound": self.results.best_objective_bound, + "wallclock_time": self.results.wallclock_time, + "expressions": self.expressions, + } + + # "best_feasible_objective", "best_objective_bound", and "wallclock_time" are all numbers, dont need subhandlers + + # subhandle "termination_condition" + results_dict["termination_condition"] = { + "value": self.results.termination_condition.value, + "name": self.results.termination_condition.name, + } + + # subhandle "solution_loader" + results_dict["solution_loader"] = {"primals": {}} + + for key, val in self.results.solution_loader.get_primals()._dict.items(): + tmp_key = key + + # handle binary vars by delving one layer in + results_dict["solution_loader"]["primals"][tmp_key] = { + "name": val[0].name, + "value": val[0].value, + "bounds": val[0].bounds, + } + + # handle binary + if val[0].is_binary(): + results_dict["solution_loader"]["primals"][tmp_key]["is_binary"] = val[ + 0 + ].is_binary() + # handle units, sometimes they dont have anything + if val[0].get_units() is not None: + results_dict["solution_loader"]["primals"][tmp_key]["units"] = ( + val[0].get_units().name + ) + else: + results_dict["solution_loader"]["primals"][tmp_key]["units"] = val[ + 0 + ].get_units() + + # renest "termination_condition" as a json-friendly dictionary + # things are either vars (which have some sort of signifier in [] brackets) or are an attribute, which dont + # the name variable will give it away + results_dict["primals_tree"] = {} + results_dict["expressions_tree"] = {} + + for key, val in self.results.solution_loader.get_primals()._dict.items(): + # split the name to figure out depth + split_name = val[0].name.split(".") + + # start at the bottom and nest accordingly + tmp_dict = { + "name": val[0].name, + "value": val[0].value, + "bounds": val[0].bounds, + } + + # handle binary + if val[0].is_binary(): + tmp_dict["is_binary"] = val[0].is_binary() + + # handle units, sometimes they dont have anything + if val[0].get_units() is not None: + tmp_dict["units"] = val[0].get_units().name + else: + tmp_dict["units"] = val[0].get_units() + + # allocate the nested dictionary + def nested_set(this_dict, key, val): + if len(key) > 1: + # check if it's a binary var and pull up one layer + if key[1] == "binary_indicator_var": + this_dict[key[0]] = val + else: + this_dict.setdefault(key[0], {}) + nested_set(this_dict[key[0]], key[1:], val) + else: + this_dict[key[0]] = val + + nested_set(results_dict["primals_tree"], split_name, tmp_dict) + + pass + + for key, val in self.expressions.items(): + # split the name to figure out depth + split_name = key.split(".") + + # start at the bottom and nest accordingly + tmp_dict = { + "value": val, + } + + # allocate the nested dictionary + def nested_set(this_dict, key, val): + if len(key) > 1: + # check if it's a binary var and pull up one layer + this_dict.setdefault(key[0], {}) + nested_set(this_dict[key[0]], key[1:], val) + else: + this_dict[key[0]] = val + + nested_set(results_dict["expressions_tree"], split_name, tmp_dict) + # split out expressions + self.expressions_tree = results_dict["expressions_tree"] + + # mint the final dictionary to save + out_dict = {"data": self.data.data, "results": results_dict} + + self.primals_tree = results_dict["primals_tree"] + + return out_dict + + def discover_level_relationships(self, dispatch_level_dict): + list_of_keys = list(dispatch_level_dict.keys()) + + relationships_dict = {} + # go through each key and split them into categories and names + # each name should have a handful of categories, which should be the same across a group of names + for this_key in list_of_keys: + # check if it has a bracketed relationship, and if it does go ahead, otherwise skip + try: + primal_category = this_key.split("[")[0] + primal_name = this_key.split("[")[1].split("]")[0] + relationships_dict.setdefault(primal_name, set()) + relationships_dict[primal_name].add(primal_category) + + except IndexError as iEx: + print( + f'[WARNING] discover_level_relationships has encountered an error: Attempted to split out {this_key}, failed with error: "{iEx}". Assigning as axuilary.' + ) + + # convert sets to frozensets to be hashable + for this_key in relationships_dict: + relationships_dict[this_key] = frozenset(relationships_dict[this_key]) + + # now go through each primal name and check for the groups who match + matching_groups_dict = {} + for this_primal_name, this_primal_set in relationships_dict.items(): + matching_groups_dict.setdefault(this_primal_set, set()) + matching_groups_dict[this_primal_set].add(this_primal_name) + + return matching_groups_dict + + def _level_relationship_dict_to_df_workhorse( + self, level_key, timeseries_dict, keys_of_interest, vars_of_interest + ): + df_data_dict = {} + units_dict = {} + # set our defaults + df_data_dict.setdefault(level_key, []) + for this_koi in keys_of_interest: + for this_voi in vars_of_interest: + df_data_dict.setdefault(f"{this_koi}_{this_voi}_value", []) + df_data_dict.setdefault(f"{this_koi}_{this_voi}_lower_bound", []) + df_data_dict.setdefault(f"{this_koi}_{this_voi}_upper_bound", []) + + # dump data into dict to read into df + for period_dict in timeseries_dict: + df_data_dict[level_key].append(period_dict["period_number"]) + for this_koi in keys_of_interest: + for this_voi in vars_of_interest: + # check if this is a variable by checking if it has a "value" + if "value" in period_dict["primals_by_name"][this_koi][this_voi]: + # if its an integer, cast it as a boolean for now + if ( + "is_binary" + in period_dict["primals_by_name"][this_koi][this_voi] + ): + if period_dict["primals_by_name"][this_koi][this_voi][ + "is_binary" + ]: + df_data_dict[f"{this_koi}_{this_voi}_value"].append( + bool( + round( + period_dict["primals_by_name"][this_koi][ + this_voi + ]["value"] + ) + ) # have to cast to int because there are floating point errors + ) + units_dict.setdefault( + f"{this_koi}_{this_voi}_value", + period_dict["primals_by_name"][this_koi][this_voi][ + "units" + ], + ) + else: + df_data_dict[f"{this_koi}_{this_voi}_value"].append( + period_dict["primals_by_name"][this_koi][this_voi][ + "value" + ] + ) + units_dict.setdefault( + f"{this_koi}_{this_voi}_value", + period_dict["primals_by_name"][this_koi][this_voi][ + "units" + ], + ) + else: + df_data_dict[f"{this_koi}_{this_voi}_value"].append( + period_dict["primals_by_name"][this_koi][this_voi][ + "value" + ] + ) + units_dict.setdefault( + f"{this_koi}_{this_voi}_value", + period_dict["primals_by_name"][this_koi][this_voi][ + "units" + ], + ) + df_data_dict[f"{this_koi}_{this_voi}_lower_bound"].append( + period_dict["primals_by_name"][this_koi][this_voi][ + "bounds" + ][0] + ) + units_dict.setdefault( + f"{this_koi}_{this_voi}_value", + period_dict["primals_by_name"][this_koi][this_voi]["units"], + ) + df_data_dict[f"{this_koi}_{this_voi}_upper_bound"].append( + period_dict["primals_by_name"][this_koi][this_voi][ + "bounds" + ][1] + ) + units_dict.setdefault( + f"{this_koi}_{this_voi}_value", + period_dict["primals_by_name"][this_koi][this_voi]["units"], + ) + + # try to make a DF, and if not just pass back an empty + try: + data_df = pd.DataFrame(df_data_dict) + # fix any Nones and make them NaNs + data_df = data_df.fillna(value=np.nan) + return data_df, units_dict + except ValueError as vEx: + print( + f"[WARNING] _level_relationship_dict_to_df_workhorse attempted to create dataframe and failed: {vEx}" + ) + return pd.DataFrame(), {} + + def _plot_workhorse_relational( + self, + level_key, + df, + keys, + vars, + parent_key_string, + pretty_title="Selected Data", + plot_bounds=False, + save_dir=".", + aspect_ratio=1, + ): + + # figure out how big the plot needs to be + gridspec_height = 2 * max(len(keys), len(vars)) + gridspec_width = 2 + fig_width_padding = 0 + fig_height_padding = 0 + max_figheight = 48 + total_periods = len(df[level_key]) + key_gridspec_div = floor( + gridspec_height / len(keys) + ) # number of gridspec heights a key plot can be + var_gridspec_div = floor( + gridspec_height / len(vars) + ) # number of gridspec heights a var plot can be + + # to make things look nice, we dont want height or width to be more than twice the other + fig_width = (total_periods * gridspec_width * 4) + fig_width_padding + fig_width = min(max_figheight, fig_width) + fig_height = (2 * gridspec_height) + fig_height_padding + if fig_width / fig_height > aspect_ratio: + fig_height = floor(fig_width / aspect_ratio) + elif fig_height / fig_width > aspect_ratio: + fig_width = floor(fig_height / aspect_ratio) + + # set up plot + fig = plt.figure( + figsize=(fig_width, fig_height), tight_layout=False + ) # (32, 16) works will for 4 plots tall and about 6 periods wide per plot + gs = fig.add_gridspec(gridspec_height, gridspec_width) + # plot out the keys of interest + ax_koi_list = [] + for ix_koi, this_koi in enumerate(keys): + ax_koi = fig.add_subplot( + gs[(ix_koi * key_gridspec_div) : ((ix_koi + 1) * key_gridspec_div), 0] + ) + ax_koi_list.append(ax_koi) + + for iy, this_voi in enumerate(vars): + ax_koi.plot( + df[level_key], + df[f"{this_koi}_{this_voi}_value"], + label=f"{this_koi}_{this_voi}", + marker="o", + ) + if plot_bounds: + ax_koi.fill_between( + df[level_key], + df[f"{this_koi}_{this_voi}_lower_bound"], + df[f"{this_koi}_{this_voi}_upper_bound"], + alpha=0.1, + ) + + ax_koi.set_ylabel("Value $[n]$") + ax_koi.xaxis.set_major_locator(MaxNLocator(integer=True)) + ax_koi.legend() + + # label axes + ax_koi_list[-1].set_xlabel(f"{level_key} $[n]$") + ax_koi_list[0].set_title(f"{pretty_title} by Type") + + # plot variables of interest + ax_voi_list = [] + # plot generations and curtailmentsagainst each outher + for ix_voi, this_voi in enumerate(vars): + ax_voi = fig.add_subplot( + gs[(ix_voi * var_gridspec_div) : ((ix_voi + 1) * var_gridspec_div), 1] + ) + ax_voi_list.append(ax_voi) + for this_koi in keys: + ax_voi.plot( + df[level_key], + df[f"{this_koi}_{this_voi}_value"], + label=f"{this_koi}_{this_voi}", + marker="o", + ) + if plot_bounds: + ax_voi.fill_between( + df[level_key], + df[f"{this_koi}_{this_voi}_lower_bound"], + df[f"{this_koi}_{this_voi}_upper_bound"], + alpha=0.1, + ) + + ax_voi.set_ylabel("Value $[n]$") + ax_voi.xaxis.set_major_locator(MaxNLocator(integer=True)) + ax_voi.legend() + + # label axes + ax_voi_list[-1].set_xlabel(f"{level_key} $[n]$") + ax_voi_list[0].set_title(f"{pretty_title} by Category") + + fig.align_labels() + fig.suptitle(f"{parent_key_string}") + fig.savefig( + f"{save_dir}{parent_key_string}_{pretty_title.replace(' ', '_')}.png" + ) + plt.close() + + def _plot_workhose_binaries( + self, + level_key, + df, + keys, + vars, + parent_key_string, + pretty_title="Selected Data", + save_dir=".", + ): + + fig = plt.figure(figsize=(32, 16), tight_layout=False) + gs = fig.add_gridspec(1, 1) # only need 1 plot for now + # if all the variables are binaries, we can assume that the vars are all binaries and the keys are all categories + total_height = len(vars) + interstate_height = 1.0 / (len(keys) + 2) + width = 1 + width_padding = 0.05 + ax_bins = fig.add_subplot(gs[:, :]) + ax_bins.set_ylim([-0.5, total_height - 0.5]) # set ylims to support bools + ax_bins.set_xlim([0.5, len(df[level_key]) + 0.5]) # set xlims to support bools + ax_bins.set_yticklabels([None] + list(vars)) + ax_bins.yaxis.set_major_locator(MaxNLocator(integer=True)) + ax_bins.xaxis.set_major_locator(MaxNLocator(integer=True)) + for axline_ix in range(total_height): + ax_bins.axhline( + axline_ix + 0.5, + color="grey", + linewidth=3, + ) # draw a seperator line between each level + for axline_ix in range(len(df[level_key])): + ax_bins.axvline( + axline_ix + 0.5, + color="grey", + linewidth=3, + linestyle="dotted", + alpha=0.5, + ) # draw a seperator line between each level + + for ix_key, this_koi in enumerate(keys): + # make a dummy line to steal the color cycler and make a single item for the legend + (line,) = ax_bins.plot( + [None], + [None], + label=f"{this_koi}", + linewidth=5, + ) + for ix_var, this_voi in enumerate(vars): + for tx, is_it_on in zip( + df[level_key], df[f"{this_koi}_{this_voi}_value"] + ): + if is_it_on: + tmp_rect = plt.Rectangle( + [ + tx - 0.5 + width_padding, + ((ix_var) + (interstate_height * (ix_key + 1))) - 0.5, + ], + width - (width_padding * 2), + interstate_height, + alpha=0.9, + edgecolor="black", + color=line.get_color(), + ) + ax_bins.add_patch(tmp_rect) + + ax_bins.set_xlabel(f"{level_key} $[n]$") + ax_bins.set_title("State Variable Time History") + ax_bins.set_ylabel("Binary State") + ax_bins.legend() + + fig.align_labels() + fig.suptitle(f"{parent_key_string}") + fig.savefig( + f"{save_dir}{parent_key_string}_{pretty_title.replace(' ', '_')}.png" + ) + plt.close() + + def _level_relationship_df_to_plot( + self, + level_key, + df, + keys, + vars, + parent_key_string, + pretty_title="Selected Data", + plot_bounds=False, + save_dir=".", + config={}, + ): + + # [HACK] hard coding the generator state order, to be fixed later + config["order_gen_state"] = ["genOff", "genShutdown", "genStartup", "genOn"] + config["order_gen_invest_state"] = [ + "genDisabled", + "genRetired", + "genExtended", + "genInstalled", + "genOperational", + ] + config["order_branch_invest_state"] = [ + "branchDisabled", + "branchRetired", + "branchExtended", + "branchInstalled", + "branchOperational", + ] + + # check if ALL the possible things to look at are binaries + all_binaries = True + for ix, this_voi in enumerate(vars): + for iy, this_koi in enumerate(keys): + if not (df[f"{this_koi}_{this_voi}_value"].dtype == "bool"): + all_binaries = False + break + if all_binaries: + + # check the config to see if we have any overrides + if "order_gen_state" in config: + # check that everything can be mapped over + matched_config_override = True + for item in vars: + if not item in config["order_gen_state"]: + matched_config_override = False + break + if matched_config_override: + vars = config["order_gen_state"] + if "order_gen_invest_state" in config: + # check that everything can be mapped over + matched_config_override = True + for item in vars: + if not item in config["order_gen_invest_state"]: + matched_config_override = False + break + if matched_config_override: + vars = config["order_gen_invest_state"] + if "order_branch_invest_state" in config: + # check that everything can be mapped over + matched_config_override = True + for item in vars: + if not item in config["order_branch_invest_state"]: + matched_config_override = False + break + if matched_config_override: + vars = config["order_branch_invest_state"] + + self._plot_workhose_binaries( + level_key, + df, + keys, + vars, + parent_key_string, + pretty_title, + save_dir, + ) + + else: + self._plot_workhorse_relational( + level_key, + df, + keys, + vars, + parent_key_string, + pretty_title, + plot_bounds, + save_dir, + ) + + def _expressions_plot_workhorse( + self, + level_key, + upper_level_dict, + parent_key_string, + save_dir="./", + plot_bounds=False, + ): + # go through a commitment period and parse out the dispatch periods + # slice out all keys pertaining to dispatchPeriod + level_period_keys = [ + this_key for this_key in upper_level_dict.keys() if (level_key in this_key) + ] + + # scan level period for keys that have values associated + keys_of_vals_of_interest = [] + for this_key in upper_level_dict[level_period_keys[0]]: + if "value" in upper_level_dict[level_period_keys[0]][this_key]: + keys_of_vals_of_interest.append(this_key) + + print(keys_of_vals_of_interest) + + # if we found things that have values, make a dictionary and then plot + # why do we need to do it by level first and then pivot the dict? Because we don't actually know what the "periods" numbers are, they might be arbitrary and in an arbitrary order + if len(keys_of_vals_of_interest) > 0: + # check all the periods and sort them out + vals_dict = {} + + for this_key in upper_level_dict: + level_period_number = int( + re.split("\[|\]", this_key.split(level_key)[1])[1] + ) + vals_dict.setdefault(level_period_number, {}) + for this_val_key in keys_of_vals_of_interest: + vals_dict[level_period_number][this_val_key] = upper_level_dict[ + this_key + ][this_val_key]["value"] + + print(vals_dict) + + # now pivot the dictionary to make a dataframe + # make a dictionary where the keys are the top layer + df_dict = {key: [] for key in keys_of_vals_of_interest} + sorted_vals_period = sorted(vals_dict) + df_dict["period_number"] = sorted_vals_period + for this_val_period in sorted_vals_period: + for this_key in keys_of_vals_of_interest: + df_dict[this_key].append(vals_dict[this_val_period][this_key]) + + print(df_dict) + expression_level_df = pd.DataFrame(df_dict) + print(expression_level_df) + + # plot the DF + # figure out how big the plot needs to be + gridspec_height = 2 * len(keys_of_vals_of_interest) + gridspec_width = 2 + fig_width_padding = 0 + fig_height_padding = 0 + max_figheight = 48 + total_periods = len(expression_level_df) + key_gridspec_div = floor( + gridspec_height / len(keys_of_vals_of_interest) + ) # number of gridspec heights a key plot can be + + # to make things look nice, we dont want height or width to be more than twice the other + fig_width = (total_periods * gridspec_width * 4) + fig_width_padding + fig_width = min(max_figheight, fig_width) + fig_height = (2 * gridspec_height) + fig_height_padding + if fig_width / fig_height > 2: + fig_height = floor(fig_width / 2) + elif fig_height / fig_width > 2: + fig_width = floor(fig_height / 2) + + # set up plot + fig = plt.figure( + figsize=(fig_width, fig_height), tight_layout=False + ) # (32, 16) works will for 4 plots tall and about 6 periods wide per plot + gs = fig.add_gridspec(gridspec_height, gridspec_width) + # plot out the keys of interest + + pretty_title = "Expression" + + ax_koi_list = [] + for ix_koi, this_koi in enumerate(keys_of_vals_of_interest): + ax_koi = fig.add_subplot( + gs[ + (ix_koi * key_gridspec_div) : ((ix_koi + 1) * key_gridspec_div), + :, + ] + ) + ax_koi_list.append(ax_koi) + + ax_koi.plot( + expression_level_df["period_number"], + expression_level_df[this_koi], + label=f"{this_koi}", + marker="o", + ) + ax_koi.set_ylabel("Value $[n]$") + ax_koi.xaxis.set_major_locator(MaxNLocator(integer=True)) + ax_koi.legend() + + # label axes + ax_koi_list[-1].set_xlabel(f"{level_key} $[n]$") + ax_koi_list[0].set_title(f"{pretty_title} by Type") + + # JSC update - " ", "_" to ' ', '_' for compilation. Not sure if this is due to a version diff or what + fig.align_labels() + fig.suptitle(f"{parent_key_string}") + fig.savefig( + f"{save_dir}{parent_key_string}_{pretty_title.replace(' ', '_')}.png" + ) + plt.close() + + def _level_plot_workhorse( + self, + level_key, + upper_level_dict, + parent_key_string, + save_dir="./", + plot_bounds=False, + ): + # go through a commitment period and parse out the dispatch periods + level_timeseries = [] + # slice out all keys pertaining to dispatchPeriod + level_period_keys = [ + this_key for this_key in upper_level_dict.keys() if (level_key in this_key) + ] + # aux_var_dict = {} + for this_key in level_period_keys: + level_period_dict = {} + # cut out which dispatch period this is + level_period_number = int( + re.split("\[|\]", this_key.split(level_key)[1])[1] + ) + # print(level_period_number) + + level_period_dict["period_number"] = level_period_number + + # pivot the dictionary to get the primals into categories + primals_by_category = {} + primals_by_name = {} + + # split key on brackets to get title + for this_primal in upper_level_dict[this_key]: + # check if it has a bracketed relationship, and if it does go ahead, otherwise skip + tmp_save_primal = upper_level_dict[this_key][this_primal] + try: + + # check if this_primal is splitable on brackets + + # based on the name of the primal, split out the "category" and the "name" + # Example: commitmentPeriod[1] + # - category: "commitmentPeriod" + # - name: "1" + primal_category = this_primal.split("[")[0] + primal_name = this_primal.split("[")[1].split("]")[0] + + # create one view that shares the categories, and one the shares the names + primals_by_category.setdefault(primal_category, {}) + primals_by_name.setdefault(primal_name, {}) + + primals_by_category[primal_category][primal_name] = tmp_save_primal + primals_by_name[primal_name][primal_category] = tmp_save_primal + + except IndexError as iEx: + print( + f"[WARNING] _level_plot_workhorse has encountered an error: Attempted to split out {this_primal} from {this_key}, failed with error {iEx}. Skipping." + ) + # aux_var_dict.setdefault(this_primal, []) + # aux_var_dict[this_primal].append(this_key) + level_period_dict["primals_by_category"] = primals_by_category + level_period_dict["primals_by_name"] = primals_by_name + level_timeseries.append(level_period_dict) + + # # clean up the aux vars based on what we found + # for aux_var_primal_name, aux_var_category_name in aux_var_dict.items(): + # print(aux_var_primal_name, aux_var_category_name) + + # sort by the dispatch period number + level_timeseries = sorted(level_timeseries, key=lambda x: x["period_number"]) + + # discover the relationships at the dispatch level + # ASSUMES that all the dispatch levels have the exact same underlying variables and relationships + level_relationships = self.discover_level_relationships( + upper_level_dict[level_period_keys[0]] + ) + # print("LEVEL RELATIONSHIPS DEBUG") + # print(level_relationships) + # print("END LEVEL RELATIONSHIPS DEBUG") + + # plot relationships + for vars_of_interest, keys_of_interest in level_relationships.items(): + + # sort the vars and keys for consistency + tmp_voi = sorted(vars_of_interest) + tmp_koi = sorted(keys_of_interest) + + # make a df for debug and also easy tabularness for plots + this_df_of_interest, this_df_units = ( + self._level_relationship_dict_to_df_workhorse( + level_key, level_timeseries, tmp_koi, tmp_voi + ) + ) + + # check if we got anything in the df + if not this_df_of_interest.empty: + + # just ram the variable names together, that'll be fine, right? + this_pretty_title = ", ".join(tmp_voi) + + # [HACK] + # if we find powerflow, plot it as a network + if "powerFlow" in tmp_voi: + self._plot_graph_workhorse( + this_df_of_interest, + "powerFlow", + parent_key_string, + units=this_df_units, + pretty_title=this_pretty_title, + save_dir=save_dir, + ) + + # [HACK] put this back one indent level when done + ## tab this back and forth to do the things + + # plot it + self._level_relationship_df_to_plot( + level_key, + this_df_of_interest, + tmp_koi, + tmp_voi, + parent_key_string, + pretty_title=this_pretty_title, + save_dir=save_dir, + plot_bounds=plot_bounds, + ) + + # JSC update - 'dc_branch' to 'branch' + def _plot_graph_workhorse( + self, + df, + value_key, + parent_key_string, + what_is_a_bus_called="branch", #'dc_branch', + units=None, + pretty_title="Selected Data", + save_dir=".", + ): + # testing networkx plots + + # preslice out data of interest + cols_of_interest = [col for col in df.columns if f"{value_key}_value" in col] + df_of_interest = df[cols_of_interest] + df_max = df_of_interest.to_numpy().max() + df_min = df_of_interest.to_numpy().min() + # assume all the units are the same, and pull the first one + units_str = "" + if units[cols_of_interest[0]] is not None: + units_str = f" [{units[cols_of_interest[0]]}]" + + # construct graph object + G = nx.Graph() + labels = {} + # add nodes + for item in self.data.data["elements"]["bus"]: + G.add_node(item) + labels[item] = item + + # do edges manually later + + # set up plot + fig = plt.figure( + figsize=(16, 8), tight_layout=False + ) # (32, 16) works will for 4 plots tall and about 6 periods wide per plot + ax_graph = fig.add_subplot() + ax_graph.grid(False) + + # # add edges + # for item in self.data.data['elements']['branch']: + # G.add_edge(self.data.data['elements']['branch'][item]['from_bus'], self.data.data['elements']['branch'][item]['to_bus']) + + # G = nx.path_graph(5) + graph_node_position_dict = nx.kamada_kawai_layout(G) + # graph_node_position_dict = nx.planar_layout(G) + # graph_node_position_dict = nx.spectral_layout(G) + nx.drawing.draw_networkx_nodes( + G, graph_node_position_dict, node_size=1000, ax=ax_graph + ) + nx.draw_networkx_labels( + G, + graph_node_position_dict, + labels, + font_size=18, + font_color="whitesmoke", + ax=ax_graph, + ) + + def draw_single_edge_flow( + item, + glyph_values_slice, + ax_graph, + cmap=cm.rainbow, + norm=Normalize(vmin=None, vmax=None), + glyph_type="custom", + ): + + def generate_flow_glyphs( + num_glyphs, + spacing=0.05, + glyph_type="triangle", + glyph_rotation=0.0, + verts=3, + ): + + flow_glyphs = [] + for this_block_ix in range(num_glyphs): + # normalizing this patch to 1 + + #### + # rectangle version + #### + if glyph_type == "rectangle": + # anchor for rectangles are set to bottom left + glyph_anchor_coord = [this_block_ix / float(num_glyphs), -0.5] + # height is y, width is x + consistent_width = 1.0 / float(num_glyphs) + # apply scaling + x_nudge = consistent_width * (spacing) + # nudge the start forward a bit (by the nudge factor) + glyph_anchor_coord[0] += x_nudge + patch_width = consistent_width - x_nudge + patch_height = 1 + flow_glyphs.append( + Rectangle(glyph_anchor_coord, patch_width, patch_height) + ) + + #### + # triangle version + #### + if glyph_type == "triangle": + # triangles need to be in the center and then given a size + glyph_anchor_coord = [ + (this_block_ix + 0.5) / float(num_glyphs), + 0, + ] + glyph_verts = 3 + glyph_radius = (1.0 / float(num_glyphs)) / 2.0 + # apply nudges + glyph_radius *= 1 - (spacing / 2.0) + flow_glyphs.append( + RegularPolygon( + glyph_anchor_coord, + glyph_verts, + radius=glyph_radius, + orientation=glyph_rotation, + ) + ) + + yscale_transform = Affine2D().scale(sx=1, sy=0.5 / glyph_radius) + # rescale y to make it fit in a 1x1 box + flow_glyphs[-1].set_transform(yscale_transform) + + #### + # n-gon version + #### + if glyph_type == "n-gon": + # triangles need to be in the center and then given a size + glyph_anchor_coord = [ + (this_block_ix + 0.5) / float(num_glyphs), + 0, + ] + glyph_verts = verts + glyph_radius = (1.0 / float(num_glyphs)) / 2.0 + # apply nudges + glyph_radius *= 1 - (spacing) + flow_glyphs.append( + RegularPolygon( + glyph_anchor_coord, + glyph_verts, + radius=glyph_radius, + orientation=glyph_rotation, + ) + ) + + yscale_transform = Affine2D().scale(sx=1, sy=0.5 / glyph_radius) + # rescale y to make it fit in a 1x1 box + flow_glyphs[-1].set_transform(yscale_transform) + + #### + # custom_flow + #### + if glyph_type == "custom": + # anchor for rectangles are set to bottom left + glyph_anchor_coord = [this_block_ix / float(num_glyphs), -0.5] + # height is y, width is x + consistent_width = 1.0 / float(num_glyphs) + # apply scaling + x_nudge = consistent_width * (spacing) + patch_width = consistent_width - x_nudge + patch_height = 1 + codes, verts = zip( + *[ + (mpath.Path.MOVETO, glyph_anchor_coord), + ( + mpath.Path.LINETO, + [ + glyph_anchor_coord[0], + glyph_anchor_coord[1] + patch_height, + ], + ), + ( + mpath.Path.LINETO, + [ + glyph_anchor_coord[0] + patch_width * 0.7, + glyph_anchor_coord[1] + patch_height, + ], + ), # go 70% of the width along the top + ( + mpath.Path.LINETO, + [ + glyph_anchor_coord[0] + patch_width, + glyph_anchor_coord[1] + patch_height * 0.5, + ], + ), # go the rest of the width and meet in the center + ( + mpath.Path.LINETO, + [ + glyph_anchor_coord[0] + patch_width * 0.7, + glyph_anchor_coord[1], + ], + ), # go back a bit and to the bottom to finish the wedge + (mpath.Path.LINETO, glyph_anchor_coord), + ] + ) # go to home + + flow_glyphs.append( + PathPatch(mpath.Path(verts, codes), ec="none"), + ) + + rotation_transofrm = Affine2D().rotate_around( + glyph_anchor_coord[0] + patch_width * 0.5, + glyph_anchor_coord[1] + patch_height * 0.5, + glyph_rotation, + ) + # rescale y to make it fit in a 1x1 box + flow_glyphs[-1].set_transform(rotation_transofrm) + + return flow_glyphs + + # make some blocks + # weights_top = (np.random.randn(4)+1)/2. + # weights_bot = (np.random.randn(4)+1)/2. + # weights_top = np.array(range(num_blocks))/(num_blocks*2) + # weights_bot = (np.array(range(num_blocks))+num_blocks)/(num_blocks*2) + weights_top = glyph_values_slice + weights_bot = glyph_values_slice + + top_flow_glyphs = generate_flow_glyphs( + len(weights_top), glyph_type=glyph_type + ) + top_facecolors = cmap(norm(weights_top)) + top_flow_collection = PatchCollection( + top_flow_glyphs, facecolors=top_facecolors, edgecolors="grey", alpha=0.5 + ) + # bot_flow_glyphs = generate_flow_glyphs(len(weights_bot), glyph_type=glyph_type, glyph_rotation=(np.pi/2.)) # for squares + bot_flow_glyphs = generate_flow_glyphs( + len(weights_bot), glyph_type=glyph_type, glyph_rotation=(np.pi) + ) # for custom + bot_flow_glyphs = reversed(bot_flow_glyphs) + bot_facecolors = cmap(norm(weights_bot)) + # bot_flow_collection = PatchCollection(bot_flow_glyphs, facecolors=bot_facecolors, edgecolors='grey', alpha=0.5) # [HACK] + + # scale and move top and bottom collections + # top_base_transform = Affine2D().scale(sx=1, sy=0.9) + Affine2D().translate(0, 0.5) #+ ax_graph.transData # [HACK] + top_base_transform = Affine2D().scale(sx=1, sy=1.0) + Affine2D().translate( + 0, 0.0 + ) # + ax_graph.transData + top_flow_collection.set_transform(top_base_transform) + bot_base_transform = Affine2D().scale(sx=1, sy=0.9) + Affine2D().translate( + 0, -0.5 + ) # + ax_graph.transData + # bot_base_transform = Affine2D().scale(sx=1, sy=0.9) + Affine2D().translate(0, -0.5) + ax_graph.transData + # bot_flow_collection.set_transform(bot_base_transform) # [HACK] + + # combine collections and move to edge between nodes + + # attempt to rotate + start_key = self.data.data["elements"][what_is_a_bus_called][item][ + "from_bus" + ] + end_key = self.data.data["elements"][what_is_a_bus_called][item]["to_bus"] + start_pos = graph_node_position_dict[start_key] + end_pos = graph_node_position_dict[end_key] + node_distance = np.linalg.norm(end_pos - start_pos) + rot_angle_rad = np.arctan2( + (end_pos[1] - start_pos[1]), (end_pos[0] - start_pos[0]) + ) + + along_edge_scale = 0.4 + away_from_edge_scale = 0.1 + # set up transformations + # stretch to the distance between target nodes + length_transform = Affine2D().scale( + sx=node_distance * along_edge_scale, sy=1 + ) + # squish + scale_transform = Affine2D().scale(sx=1, sy=away_from_edge_scale) + # rotate + rot_transform = Affine2D().rotate_deg(np.rad2deg(rot_angle_rad)) + # translate to the node start, then push it along the edge until it's apprximately centered and scaled nicely + translate_transform = Affine2D().translate( + start_pos[0] + + ( + np.cos(rot_angle_rad) * node_distance * 0.5 * (1 - along_edge_scale) + ), + start_pos[1] + + ( + np.sin(rot_angle_rad) * node_distance * 0.5 * (1 - along_edge_scale) + ), + ) + t2 = ( + length_transform + + scale_transform + + rot_transform + + translate_transform + + ax_graph.transData + ) + + top_flow_collection.set_transform(top_flow_collection.get_transform() + t2) + # bot_flow_collection.set_transform(bot_flow_collection.get_transform() + t2) # [HACK] + + # add collection + ax_graph.add_collection(top_flow_collection) + # ax_graph.add_collection(bot_flow_collection) # [HACK] + + # add edges + # define edge colorbar + cmap = cm.rainbow + normalize = Normalize(vmin=df_min, vmax=df_max) + cmappable = cm.ScalarMappable(norm=normalize, cmap=cmap) + + for item in self.data.data["elements"][what_is_a_bus_called]: + + # grab the keys we care about + start_key = self.data.data["elements"][what_is_a_bus_called][item][ + "from_bus" + ] + end_key = self.data.data["elements"][what_is_a_bus_called][item]["to_bus"] + start_pos = graph_node_position_dict[start_key] + end_pos = graph_node_position_dict[end_key] + # edge_key = f"branch_{start_key}_{end_key}_{value_key}_value" + # alt_edge_key = f"branch_{end_key}_{start_key}_{value_key}_value" + edge_key = f"{item}_{value_key}_value" + alt_edge_key = f"{item}_{value_key}_value" + + # @KyleSkolfield is there a reason to not do this? + branch_name_edge_key = item + "_powerFlow_value" + + # kind = 'triangle' + # kind = 'rectangle' + kind = "custom" + glyph_values_slice = df[branch_name_edge_key].values + + try: + glyph_values_slice = df[edge_key].values + except KeyError as kex: + print(f"Attempted to slice DF in network using {edge_key}, failed.") + + draw_single_edge_flow( + item, + glyph_values_slice, + ax_graph, + cmap=cmap, + norm=normalize, + glyph_type=kind, + ) + + # forward arrow + ax_graph.arrow( + start_pos[0], + start_pos[1], + (end_pos[0] - start_pos[0]), + (end_pos[1] - start_pos[1]), + color="black", + ) + # backward arrow + ax_graph.arrow( + end_pos[0], + end_pos[1], + (start_pos[0] - end_pos[0]), + (start_pos[1] - end_pos[1]), + color="black", + ) + + # insert colorbar + fig.colorbar(cmappable, ax=ax_graph, label=f"{value_key}{units_str}") + # make some titles + fig.suptitle(f"{parent_key_string}_{value_key}") + + # JSC update - " ", "_" to ' ', '_' for compilation. Not sure if this is due to a version diff or what + # save + fig.savefig( + f"{save_dir}{parent_key_string}_{pretty_title.replace(' ', '_')}_graph.png" + ) + pass + + def plot_levels(self, save_dir="."): + + # self._plot_graph_workhorse() + + # plot or represent primals trees + for this_root_level_key in self.primals_tree: + if "investmentStage" in this_root_level_key: + # run the toplevel keys + parent_key_string = f"{this_root_level_key}" + self._level_plot_workhorse( + "investmentStage", self.primals_tree, this_root_level_key, save_dir + ) + + # run the representative period subkeys + investment_level_cut = self.primals_tree[this_root_level_key] + parent_key_string = f"{this_root_level_key}" + self._level_plot_workhorse( + "representativePeriod", + investment_level_cut, + parent_key_string, + save_dir, + ) + + for this_inv_level_key in self.primals_tree[this_root_level_key].keys(): + if "representativePeriod" in this_inv_level_key: + representative_level_cut = self.primals_tree[ + this_root_level_key + ][this_inv_level_key] + parent_key_string = ( + f"{this_root_level_key}_{this_inv_level_key}" + ) + self._level_plot_workhorse( + "commitmentPeriod", + representative_level_cut, + parent_key_string, + save_dir, + ) + + for this_rep_level_key in self.primals_tree[ + this_root_level_key + ][this_inv_level_key].keys(): + if "commitmentPeriod" in this_rep_level_key: + commitment_level_cut = self.primals_tree[ + this_root_level_key + ][this_inv_level_key][this_rep_level_key] + + parent_key_string = f"{this_root_level_key}_{this_inv_level_key}_{this_rep_level_key}" + + self._level_plot_workhorse( + "dispatchPeriod", + commitment_level_cut, + parent_key_string, + save_dir, + ) + + # # plot or represent expressions + # self._expressions_plot_workhorse( + # "investmentStage", self.expressions_tree, 'investmentStage', save_dir + # ) + # for this_root_level_key in self.expressions_tree: + # if "investmentStage" in this_root_level_key: + # investment_level_cut = self.expressions_tree[this_root_level_key] + # parent_key_string = f"{this_root_level_key}" + # self._expressions_plot_workhorse( + # "representativePeriod", investment_level_cut, parent_key_string, save_dir + # ) + + # for this_inv_level_key in self.expressions_tree[this_root_level_key].keys(): + # if "representativePeriod" in this_inv_level_key: + # representative_level_cut = self.expressions_tree[this_root_level_key][this_inv_level_key] + # parent_key_string = f"{this_root_level_key}_{this_inv_level_key}" + # self._expressions_plot_workhorse( + # "commitmentPeriod", representative_level_cut, parent_key_string, save_dir + # ) + + # for this_rep_level_key in self.expressions_tree[ + # this_root_level_key + # ][this_inv_level_key].keys(): + # if "commitmentPeriod" in this_rep_level_key: + # commitment_level_cut = self.expressions_tree[ + # this_root_level_key + # ][this_inv_level_key][this_rep_level_key] + + # parent_key_string = f"{this_root_level_key}_{this_inv_level_key}_{this_rep_level_key}" + + # self._expressions_plot_workhorse( + # "dispatchPeriod",commitment_level_cut, parent_key_string, save_dir + # ) + # pass diff --git a/gtep/driver_idaes_viz_reliability.py b/gtep/driver_idaes_viz_reliability.py new file mode 100644 index 00000000..ce71ab90 --- /dev/null +++ b/gtep/driver_idaes_viz_reliability.py @@ -0,0 +1,64 @@ +from gtep.gtep_model import ExpansionPlanningModel +from gtep.gtep_data import ExpansionPlanningData +from gtep.gtep_solution import ExpansionPlanningSolution +from pyomo.core import TransformationFactory +from pyomo.contrib.appsi.solvers.highs import Highs +from pyomo.contrib.appsi.solvers.gurobi import Gurobi +import gtep.validation + + +data_path = "./gtep/data/SanDiego" +data_object = ExpansionPlanningData() +data_object.load_prescient(data_path) + +sol_object = ExpansionPlanningSolution() +sol_object.import_data_object(data_object) +sol_object.read_json("C:/Users/jkskolf/idaes-gtep/gtep_pre_reliability_solution.json") + + +data_out_path = "./Sim Eng Viz/PreReliability/Prescient/data" + +gtep.validation.populate_generators(data_path, sol_object, data_out_path) +gtep.validation.populate_transmission(data_path, sol_object, data_out_path) +gtep.validation.filter_pointers(data_path, data_out_path) +gtep.validation.clone_timeseries(data_path, data_out_path) + + +# sol_object.import_data_object(data_object) + +# sol_object.plot_levels(save_dir="./Sim Eng Viz/Hourly/GTEP/plots/") + +from prescient.simulator import Prescient + +# set some options +prescient_options = { + "data_path": data_out_path, + "input_format": "rts-gmlc", + "simulate_out_of_sample": False, + "run_sced_with_persistent_forecast_errors": False, + "output_directory": "./Sim Eng Viz/PreReliability/Prescient/results", + "start_date": "07-26-2020", + "num_days": 21, + "sced_horizon": 24, + "ruc_mipgap": 0.01, + "reserve_factor": 0, + "deterministic_ruc_solver": "gurobi_persistent", + "deterministic_ruc_solver_options": { + "feas": "off", + "DivingF": "on", + }, + "sced_solver": "gurobi", + "sced_frequency_minutes": 60, + "ruc_horizon": 48, + "compute_market_settlements": True, + "monitor_all_contingencies": False, + "output_solver_logs": False, + "price_threshold": 1000, + "contingency_price_threshold": 100, + "reserve_price_threshold": 5, + "ruc_network_type": "btheta", + "sced_network_type": "btheta", + "enforce_sced_shutdown_ramprate": False, +} +# run the simulator +Prescient().simulate(**prescient_options) diff --git a/gtep/gtep_data.py b/gtep/gtep_data.py index c21b6e66..7f540422 100644 --- a/gtep/gtep_data.py +++ b/gtep/gtep_data.py @@ -102,6 +102,7 @@ def load_default_data_settings(self): self.md.data["elements"]["generator"][gen]["emissions_factor"] = 1 self.md.data["elements"]["generator"][gen]["start_fuel"] = 1 self.md.data["elements"]["generator"][gen]["investment_cost"] = 235164 + self.md.data["elements"]["generator"][gen]["failure_rate"] = 0.1 for branch in self.md.data["elements"]["branch"]: self.md.data["elements"]["branch"][branch]["loss_rate"] = 0 self.md.data["elements"]["branch"][branch]["distance"] = 1 diff --git a/gtep/gtep_model.py b/gtep/gtep_model.py index e0fc66e8..28f80a76 100644 --- a/gtep/gtep_model.py +++ b/gtep/gtep_model.py @@ -20,7 +20,7 @@ from math import ceil -from config_options import _get_model_config, _add_common_configs +from gtep.config_options import _get_model_config, _add_common_configs # Define what a USD is for pyomo units purposes @@ -1422,7 +1422,7 @@ def representative_period_rule(b, representative_period): i_s = b.parent_block() b.currentPeriod = representative_period - if m.config["include_commitment"] or m.config["include_redispatch"]: + if m.config.get("include_commitment") or m.config.get("include_redispatch"): b.commitmentPeriods = RangeSet(m.numCommitmentPeriods[representative_period]) b.commitmentPeriod = Block(b.commitmentPeriods, rule=commitment_period_rule) diff --git a/gtep/gtep_solution.py b/gtep/gtep_solution.py index 7fcf67c3..da7db8c8 100644 --- a/gtep/gtep_solution.py +++ b/gtep/gtep_solution.py @@ -40,12 +40,12 @@ def load_from_file(self): pass def load_from_model(self, gtep_model): - if type(gtep_model) is not ExpansionPlanningModel: - logger.warning( - f"Solutions must be loaded from ExpansionPlanningModel objects, not %s" - % type(gtep_model) - ) - raise ValueError + # if type(gtep_model) is not ExpansionPlanningModel: + # logger.warning( + # f"Solutions must be loaded from ExpansionPlanningModel objects, not %s" + # % type(gtep_model) + # ) + # raise ValueError if gtep_model.results is None: raise ValueError( "ExpansionPlanningSolution objects loaded from model must have a results component."