for this week

Author: Kyle Skolfield

Diff for: gtep/driver_idaes_viz.py

@@ -0,0 +1,85 @@
+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/5bus_jsc_flat"
+data_object = ExpansionPlanningData()
+data_object.load_prescient(data_path)
+
+
+mod_object = ExpansionPlanningModel(
+    stages=2,
+    data=data_object.representative_data,
+    num_reps=4,
+    len_reps=1,
+    num_commit=6,
+    num_dispatch=4,
+)
+
+
+mod_object.config["include_investment"] = True
+mod_object.config["dispatch_randomization"] = False
+mod_object.config["scale_loads"] = False
+mod_object.create_model()
+
+
+TransformationFactory("gdp.bound_pretransformation").apply_to(mod_object.model)
+TransformationFactory("gdp.bigm").apply_to(mod_object.model)
+opt = Gurobi()
+mod_object.results = opt.solve(mod_object.model)
+
+sol_object = ExpansionPlanningSolution()
+sol_object.load_from_model(mod_object)
+sol_object.dump_json("./Sim Eng Viz/4Hourly/GTEP/idaes_solution.json")
+
+data_out_path = "./Sim Eng Viz/4Hourly/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/4Hourly/Prescient/results/04-01-2020",
+    "start_date": "03-26-2020",
+    "num_days": 7,
+    "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 for: gtep/gtep_data.py

@@ -71,16 +71,16 @@ def load_prescient(self, data_path, options_dict=None):
         # default here allows up to 24 hours for periods
         time_keys = self.md.data["system"]["time_keys"]
         key_idx = time_keys.index("2020-01-01 00:00")
-        time_key_set = time_keys[key_idx : key_idx + 24]
+        time_key_set = time_keys[key_idx : key_idx + 24 : 4]
         data_list.append(self.md.clone_at_time_keys(time_key_set))
         key_idx = time_keys.index("2020-04-01 00:00")
-        time_key_set = time_keys[key_idx : key_idx + 24]
+        time_key_set = time_keys[key_idx : key_idx + 24 : 4]
         data_list.append(self.md.clone_at_time_keys(time_key_set))
         key_idx = time_keys.index("2020-07-01 00:00")
-        time_key_set = time_keys[key_idx : key_idx + 24]
+        time_key_set = time_keys[key_idx : key_idx + 24 : 4]
         data_list.append(self.md.clone_at_time_keys(time_key_set))
         key_idx = time_keys.index("2020-10-01 00:00")
-        time_key_set = time_keys[key_idx : key_idx + 24]
+        time_key_set = time_keys[key_idx : key_idx + 24 : 4]
         data_list.append(self.md.clone_at_time_keys(time_key_set))
 
         self.representative_data = data_list
@@ -90,8 +90,8 @@ def load_default_data_settings(self):
         """Fills in necessary but unspecified data information."""
         for gen in self.md.data["elements"]["generator"]:
             self.md.data["elements"]["generator"][gen]["lifetime"] = 3
-            self.md.data["elements"]["generator"][gen]["spinning_reserve_frac"] = 0.1
-            self.md.data["elements"]["generator"][gen]["quickstart_reserve_frac"] = 0.1
+            self.md.data["elements"]["generator"][gen]["spinning_reserve_frac"] = 0
+            self.md.data["elements"]["generator"][gen]["quickstart_reserve_frac"] = 0
             self.md.data["elements"]["generator"][gen]["capital_multiplier"] = 1
             self.md.data["elements"]["generator"][gen]["extension_multiplier"] = 1
             self.md.data["elements"]["generator"][gen]["max_operating_reserve"] = 1
@@ -101,9 +101,9 @@ def load_default_data_settings(self):
             self.md.data["elements"]["generator"][gen]["ramp_down_rate"] = 0.1
             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]["investment_cost"] = 235
         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
-        self.md.data["system"]["min_operating_reserve"] = 0.1
-        self.md.data["system"]["min_spinning_reserve"] = 0.1
+        self.md.data["system"]["min_operating_reserve"] = 0
+        self.md.data["system"]["min_spinning_reserve"] = 0

Diff for: gtep/gtep_model.py

@@ -102,7 +102,9 @@ def create_model(self):
             # 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])
+            # TODO god what have I done
+            # m.md = scale_ModelData_to_pu(self.data[0])
+            m.md = self.data[0]
         else:
             # If self.data is an Egret model data object, representative periods will just copy it unchanged
             m.data_list = None
@@ -501,6 +503,8 @@ def renewable_generation_limits(b, renewableGen):
 
     # Define bounds on renewable generator curtailment
     def curtailment_limits(b, renewableGen):
+        if "HYDRO" in renewableGen:
+            return (0, 0)
         return (0, m.renewableCapacity[renewableGen])
 
     b.renewableCurtailment = Var(
@@ -712,8 +716,9 @@ def add_dispatch_constraints(b, disp_per):
     r_p = c_p.parent_block()
     i_p = r_p.parent_block()
 
-    for key in m.loads.keys():
-        m.loads[key] *= max(0, m.rng.normal(0.5, 0.2))
+    if m.config["dispatch_randomization"]:
+        for key in m.loads.keys():
+            m.loads[key] *= max(0, m.rng.normal(0.5, 0.2))
 
     # Energy balance constraint
     @b.Constraint(m.buses)
@@ -1641,6 +1646,7 @@ def model_data_references(m):
     }
 
     # Demand at each bus
+    ## FIXME this just builds a full 8760 of loads which is not really what we want
     m.loads = {
         m.md.data["elements"]["load"][load_n]["bus"]: m.md.data["elements"]["load"][
             load_n
@@ -1693,7 +1699,7 @@ def model_data_references(m):
     ## 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 / u.hr)
+    m.fixedOperatingCost = Param(m.generators, default=0, units=u.USD / u.hr)
     m.deficitPenalty = Param(m.stages, default=1, units=u.USD / (u.MW * u.hr))
 
     # Amount of fuel required to be consumed for startup process for each generator
@@ -1722,7 +1728,7 @@ def model_data_references(m):
     # NOTE: what should this be valued at?  This being both curtailment and load shed.
     # TODO: update valuations
     m.curtailmentCost = Param(
-        initialize=2 * max(value(item) for item in m.fuelCost.values()),
+        initialize=10000 * max(value(item) for item in m.fuelCost.values()),
         units=u.USD / (u.MW * u.hr),
     )
     m.loadShedCost = Param(

Diff for: gtep/gtep_solution.py

@@ -193,8 +193,9 @@ def nested_set(this_dict, key, val):
         self.expressions_tree = results_dict["expressions_tree"]
 
         # mint the final dictionary to save
-        out_dict = {"data": self.data.data, "results": results_dict}
-
+        # out_dict = {"data": self.data.data, "results": results_dict}
+        out_dict = {"data": self.data[1].data, "results": results_dict}
+        print(self.data)
         self.primals_tree = results_dict["primals_tree"]
 
         return out_dict

Diff for: gtep/idaes_parse.py

@@ -0,0 +1,121 @@
+from gtep.gtep_solution import ExpansionPlanningSolution
+import pandas as pd
+import json
+import re
+import matplotlib as mpl
+from matplotlib import colormaps
+import seaborn as sns
+
+
+with open("./Sim Eng Viz/4Hourly/GTEP/idaes_solution.json") as fil:
+    solutions = json.load(fil)
+
+target_investment_stage = "2"
+target_representative_period = "2"
+
+filtered_solns = solutions["results"]["primals_tree"]["investmentStage[2]"][
+    "representativePeriod[2]"
+]
+
+filtered_solns.pop("_logical_to_disjunctive")
+
+
+def name_filter(gen_name):
+    return re.search(r"\[.*\]", gen_name).group(0)[1:-1]
+
+
+thermal_gen = {}
+curtailment = {}
+load_shed = {}
+for key in filtered_solns.keys():
+    for k, v in filtered_solns[key]["dispatchPeriod[1]"].items():
+        if "Gen" in k:
+            k = name_filter(k)
+            if thermal_gen.get(k):
+                thermal_gen[k].append(v["value"])
+            else:
+                thermal_gen[k] = [v["value"]]
+        if "Curt" in k:
+            if curtailment.get(k):
+                curtailment[k].append(v["value"])
+            else:
+                curtailment[k] = [v["value"]]
+        if "Shed" in k:
+            if load_shed.get(k):
+                load_shed[k].append(v["value"])
+            else:
+                load_shed[k] = [v["value"]]
+    # only for the four hourly...
+    if filtered_solns[key].get("dispatchPeriod[2]"):
+        for k, v in filtered_solns[key]["dispatchPeriod[2]"].items():
+            if "Gen" in k:
+                k = name_filter(k)
+                if thermal_gen.get(k):
+                    thermal_gen[k].append(v["value"])
+                else:
+                    thermal_gen[k] = [v["value"]]
+        for k, v in filtered_solns[key]["dispatchPeriod[3]"].items():
+            if "Gen" in k:
+                k = name_filter(k)
+                if thermal_gen.get(k):
+                    thermal_gen[k].append(v["value"])
+                else:
+                    thermal_gen[k] = [v["value"]]
+        for k, v in filtered_solns[key]["dispatchPeriod[4]"].items():
+            if "Gen" in k:
+                k = name_filter(k)
+                if thermal_gen.get(k):
+                    thermal_gen[k].append(v["value"])
+                else:
+                    thermal_gen[k] = [v["value"]]
+
+
+print(thermal_gen)
+print(curtailment)
+print(load_shed)
+
+nonzero_gennames = ["10_STEAM", "3_CT", "4_STEAM"]
+for gen in nonzero_gennames:
+    thermal_gen.pop(gen)
+
+gen_df = pd.DataFrame.from_dict(dict(sorted(thermal_gen.items())))
+print(gen_df)
+prescient_gen_df = pd.read_csv(
+    "./Sim Eng Viz/4Hourly/Prescient/results/generation_out.csv",
+    names=["gen", "val"],
+)
+
+print(prescient_gen_df)
+for label, content in prescient_gen_df.items():
+    print(f"label: {label}")
+    print(f"content: {content}", sep="\n")
+
+prescient_gen = {}
+for row in prescient_gen_df.iterrows():
+    gen = row[-1]["gen"]
+    val = row[-1]["val"]
+    if prescient_gen.get(gen):
+        prescient_gen[gen].append(val)
+    else:
+        prescient_gen[gen] = [val]
+
+for gen in nonzero_gennames:
+    prescient_gen.pop(gen)
+prescient_sorted_gen_df = pd.DataFrame.from_dict(dict(sorted(prescient_gen.items())))
+
+my_cmap = sns.color_palette("muted")
+
+
+ax1 = gen_df.plot.bar(stacked=True, color=my_cmap)
+ax1.set_title("GTEP, 4-Hourly Commitment (04/01)")
+ax1.set_xlabel("Hour")
+ax1.set_ylabel("Generation (MW)")
+ax1.figure.savefig("4blah.png")
+
+ax2 = prescient_sorted_gen_df.plot.bar(stacked=True, color=my_cmap)
+ax2.set_title("Prescient (04/01)")
+ax2.set_xlabel("Hour")
+ax2.set_ylabel("Generation (MW)")
+ax2.figure.savefig("4blah2.png")
+
+print(gen_df)

Diff for: gtep/validation.py

@@ -29,16 +29,17 @@ def gen_name_filter(gen_name):
         )
 
     solution_dict = sol_object._to_dict()["results"]["primals_tree"]
-    end_investment_stage = list(solution_dict.keys())[0]
+    end_investment_stage = list(solution_dict.keys())[-2]
     end_investment_solution_dict = {
         k: v["value"]
         for k, v in solution_dict[end_investment_stage].items()
         if gen_name_filter(k) and v["value"] > 0.5
     }
-    end_investment_gens = [
+    end_investment_thermal_gens = [
         re.search(r"\[.*\]", k).group(0)[1:-1]
         for k in end_investment_solution_dict.keys()
     ]
+    end_investment_gens = end_investment_thermal_gens
 
     # for renewable:
     # total capacity should be installed + operational + extended values
@@ -71,6 +72,18 @@ def renewable_name_filter(gen_name):
         os.makedirs(data_output_path)
     output_df.to_csv(os.path.join(data_output_path, "gen.csv"), index=False)
 
+    # We also need to update initial status values for thermal gens that have been invested
+    # By default we will set these to (-24, 0)
+    # initial_df = pd.read_csv(os.path.join(data_input_path, "initial_status.csv"))
+    # for gen in end_investment_thermal_gens:
+    #     print(gen)
+    #     if gen not in initial_df.columns:
+    #         # initial_df.insert(-1, gen, [-24, 0])
+    #         initial_df[gen] = [-24, 0]
+    #         print(initial_df)
+
+    # initial_df.to_csv(os.path.join(data_output_path, "initial_status.csv"), index=False)
+
 
 def populate_transmission(data_input_path, sol_object, data_output_path):
     # load existing and candidate generators from initial prescient data
@@ -115,6 +128,8 @@ def filter_pointers(data_input_path, data_output_path):
     # keep generators that exist at the final investment stage and remove the rest
     # keep all non-generator timeseries pointers
     matching_gen_list = [gen for gen in output_generators_df["GEN UID"]]
+
+    ## FIXME: this won't work when completely NEW renewable sites are an option
     output_df = input_pointers_df[
         input_pointers_df["Object"].isin(matching_gen_list)
         | input_pointers_df["Category"]
@@ -135,9 +150,13 @@ def clone_timeseries(data_input_path, data_output_path):
     file_list.remove("timeseries_pointers.csv")
     file_list.remove("gen.csv")
     file_list.remove("branch.csv")
+    file_list.remove("initial_status.csv")
 
     # @jkskolf, I don't think I like this ...
     for fname in file_list:
         from_file = os.path.join(data_input_path, fname)
         to_file = os.path.join(data_output_path, fname)
         shutil.copy(from_file, to_file)
+
+
+# TODO: synchronize load scaling