Merge pull request #219 from Breakthrough-Energy/daniel/hifld_transmi…

…ssion_params

feat: add functions to estimate branch parameters (impedance and rating)

prereise/gather/griddata/hifld/const.py

@@ -118,3 +118,67 @@
     "345": 345,
     "500": 500,
 }
+
+line_reactance_per_mile = {  # per-unit
+    69: 0.0096,
+    100: 0.0063,
+    115: 0.0039,
+    138: 0.0026,
+    161: 0.0021,
+    230: 0.0011,
+    345: 0.0005,
+    500: 0.0002,
+    765: 0.0001,
+}
+line_rating_short = {  # MVA
+    69: 86,
+    100: 181,
+    115: 239,
+    138: 382,
+    161: 446,
+    230: 797,
+    345: 1793,
+    500: 2598,
+    765: 5300,
+}
+line_rating_short_threshold = 50  # miles
+line_rating_surge_impedance_loading = {  # MVA
+    69: 13,
+    100: 30,
+    115: 35,
+    138: 50,
+    161: 69,
+    230: 140,
+    345: 375,
+    500: 1000,
+    765: 2250,
+}
+line_rating_surge_impedance_coefficient = 43.261
+line_rating_surge_impedance_exponent = -0.6678
+
+transformer_reactance = {  # per-unit
+    (69, 115): 0.14242,
+    (69, 138): 0.10895,
+    (69, 161): 0.14943,
+    (69, 230): 0.09538,
+    (69, 345): 0.08896,
+    (115, 161): 0.04516,
+    (115, 230): 0.04299,
+    (115, 345): 0.04020,
+    (115, 500): 0.06182,
+    (138, 161): 0.02818,
+    (138, 230): 0.03679,
+    (138, 345): 0.03889,
+    (138, 500): 0.03279,
+    (138, 765): 0.02284,
+    (161, 230): 0.06539,
+    (161, 345): 0.03293,
+    (161, 500): 0.06978,
+    (230, 345): 0.02085,
+    (230, 500): 0.01846,
+    (230, 765): 0.01616,
+    (345, 500): 0.01974,
+    (345, 765): 0.01625,
+    (500, 765): 0.00436,
+}
+transformer_rating = 800  # MVA

prereise/gather/griddata/hifld/data_process/tests/test_transmission.py

@@ -2,10 +2,13 @@
 import pytest
 from pandas.testing import assert_frame_equal, assert_series_equal
 
+from prereise.gather.griddata.hifld import const
 from prereise.gather.griddata.hifld.data_process.transmission import (
     augment_line_voltages,
     create_buses,
     create_transformers,
+    estimate_branch_impedance,
+    estimate_branch_rating,
 )
 
 
@@ -116,3 +119,68 @@ def test_create_transformers():
     )
     transformers = create_transformers(bus)
     assert_frame_equal(transformers, expected_transformers)
+
+
+def test_estimate_branch_impedance_lines():
+    branch = pd.DataFrame(
+        {"VOLTAGE": [69, 70, 345], "type": ["Line"] * 3, "length": [10, 15, 20]}
+    )
+    x = estimate_branch_impedance(branch.iloc[0], pd.Series())
+    assert x == const.line_reactance_per_mile[69] * 10
+    x = estimate_branch_impedance(branch.iloc[1], pd.Series())
+    assert x == const.line_reactance_per_mile[69] * 15
+    x = estimate_branch_impedance(branch.iloc[2], pd.Series())
+    assert x == const.line_reactance_per_mile[345] * 20
+
+
+def test_estimate_branch_impedance_transformers():
+    transformers = pd.DataFrame(
+        {"from_bus_id": [0, 1, 2], "to_bus_id": [1, 2, 3], "type": ["Transformer"] * 3}
+    )
+    bus_voltages = pd.Series([69, 230, 350, 500])
+    x = estimate_branch_impedance(transformers.iloc[0], bus_voltages)
+    assert x == const.transformer_reactance[(69, 230)]
+    x = estimate_branch_impedance(transformers.iloc[1], bus_voltages)
+    assert x == const.transformer_reactance[(230, 345)]
+    x = estimate_branch_impedance(transformers.iloc[2], bus_voltages)
+    assert x == const.transformer_reactance[(345, 500)]
+
+
+def test_estimate_branch_rating_lines():
+    branch = pd.DataFrame(
+        {
+            "VOLTAGE": [69, 140, 345, 499],
+            "type": ["Line"] * 4,
+            "length": [10, 50, 100, 150],
+        }
+    )
+    rating = estimate_branch_rating(branch.iloc[0], pd.Series())
+    assert rating == const.line_rating_short[69]
+    rating = estimate_branch_rating(branch.iloc[1], pd.Series())
+    assert rating == const.line_rating_short[138]
+    rating = estimate_branch_rating(branch.iloc[2], pd.Series())
+    assert rating == (
+        const.line_rating_surge_impedance_loading[345]
+        * const.line_rating_surge_impedance_coefficient
+        * 100 ** const.line_rating_surge_impedance_exponent
+    )
+    rating = estimate_branch_rating(branch.iloc[3], pd.Series())
+    assert rating == (
+        const.line_rating_surge_impedance_loading[500]
+        * const.line_rating_surge_impedance_coefficient
+        * 150 ** const.line_rating_surge_impedance_exponent
+    )
+
+
+def test_estimate_branch_rating_transformers():
+    transformers = pd.DataFrame(
+        {"from_bus_id": [0, 1, 2], "to_bus_id": [1, 2, 3], "type": ["Transformer"] * 3}
+    )
+    bus_voltages = pd.Series([69, 230, 350, 500])
+
+    rating = estimate_branch_rating(transformers.iloc[0], bus_voltages)
+    assert rating == const.transformer_rating
+    rating = estimate_branch_rating(transformers.iloc[1], bus_voltages)
+    assert rating == const.transformer_rating * 3
+    rating = estimate_branch_rating(transformers.iloc[2], bus_voltages)
+    assert rating == const.transformer_rating * 4

prereise/gather/griddata/hifld/data_process/transmission.py

@@ -1,4 +1,6 @@
+import hashlib
 import os
+import pickle
 
 import networkx as nx
 import numpy as np
@@ -13,6 +15,9 @@
     get_hifld_electric_substations,
     get_zone,
 )
+from prereise.gather.griddata.hifld.data_process.topology import (
+    connect_islands_with_minimum_cost,
+)
 
 
 def check_for_location_conflicts(substations):
@@ -319,6 +324,7 @@ def map_via_neighbor_voltages(lines, neighbors, func, method_name):
         while True:
             missing = lines.query("VOLTAGE.isnull()")
             if len(missing) == 0:
+                print(f"No more missing voltages remain after neighbor {method_name}")
                 break
             found_voltages = missing.apply(
                 lambda x: {
@@ -433,6 +439,118 @@ def create_transformers(bus):
     return pd.DataFrame(bus_pairs, columns=["from_bus_id", "to_bus_id"])
 
 
+def estimate_branch_impedance(branch, bus_voltages):
+    """Estimate branch impedance using transformer voltages or line voltage and length.
+
+    :param pandas.Series branch: data for a single branch (line or transformer). All
+        branches require 'type' attributes, lines require 'VOLTAGE' and 'length',
+        transformers require 'from_bus_id' and 'to_bus_id'.
+    :param pandas.Series bus_voltages: mapping of buses to voltages.
+    :return: (*float*) -- impedance for that branch (per-unit).
+    """
+
+    def _euclidian(a, b):
+        return ((a[0] - b[0]) ** 2 + (a[1] - b[1]) ** 2) ** (1 / 2)
+
+    if branch.loc["type"] == "Transformer":
+        voltage_tuple = sorted(bus_voltages.loc[[branch.from_bus_id, branch.to_bus_id]])
+        reactance_lookup = pd.Series(const.transformer_reactance)
+        # Find the 'closest' voltage pair via Euclidian distance
+        closest_voltage_tuple = (
+            reactance_lookup.index.to_series()
+            .map(lambda x: _euclidian(x, voltage_tuple))
+            .idxmin()
+        )
+        return const.transformer_reactance[closest_voltage_tuple]
+    elif branch.loc["type"] == "Line":
+        # Calculate line length via sum of distances between adjacent coordinate pairs
+        reactance_lookup = pd.Series(const.line_reactance_per_mile)
+        closest_voltage_reactance_per_mile = reactance_lookup.iloc[
+            reactance_lookup.index.get_loc(branch.loc["VOLTAGE"], method="nearest")
+        ]
+        return branch.loc["length"] * closest_voltage_reactance_per_mile
+    else:
+        raise ValueError(f"{branch.loc['type']} not a valid branch type")
+
+
+def calculate_branch_mileage(branch):
+    """Estimate distance of a line.
+
+    :param pandas.Series branch: data for a single line.
+    :return: (*float*) -- distance (in miles) for that line.
+    """
+    coordinates = branch.loc["COORDINATES"]
+    return sum([haversine(a, b) for a, b in zip(coordinates[:-1], coordinates[1:])])
+
+
+def estimate_branch_rating(branch, bus_voltages):
+    """Estimate branch rating using line voltage or constant value for transformers.
+
+    :param pandas.Series branch: data for a single branch (line or transformer). All
+        branches require 'type' attributes, lines require 'VOLTAGE' and 'length'.
+    :param pandas.Series bus_voltages: mapping of buses to voltages.
+    :return: (*float*) -- rating for that branch (MW).
+    :raises ValueError: if branch 'type' attribute not recognized.
+    """
+    if branch.loc["type"] == "Line":
+        if branch.loc["length"] <= const.line_rating_short_threshold:
+            rating_lookup = pd.Series(const.line_rating_short)
+            closest_rating = rating_lookup.iloc[
+                rating_lookup.index.get_loc(branch.loc["VOLTAGE"], method="nearest")
+            ]
+            return closest_rating
+        else:
+            sil_lookup = pd.Series(const.line_rating_surge_impedance_loading)
+            closest_sil = sil_lookup.iloc[
+                sil_lookup.index.get_loc(branch.loc["VOLTAGE"], method="nearest")
+            ]
+            return (
+                closest_sil
+                * const.line_rating_surge_impedance_coefficient
+                * branch.loc["length"] ** const.line_rating_surge_impedance_exponent
+            )
+    elif branch.loc["type"] == "Transformer":
+        rating = const.transformer_rating
+        max_voltage = bus_voltages.loc[[branch.from_bus_id, branch.to_bus_id]].max()
+        line_capacities = pd.Series(const.line_rating_short)
+        closest_voltage_rating = line_capacities.iloc[
+            line_capacities.index.get_loc(max_voltage, method="nearest")
+        ]
+        num_addl_transformers = int(closest_voltage_rating / rating)
+        return rating * (1 + num_addl_transformers)
+    raise ValueError(f"{branch.loc['type']} not a valid branch type")
+
+
+def get_mst_edges(lines, substations):
+    """Get the set of lines which connected the connected components of the lines graph,
+    either from a cache or by generating from scratch and caching.
+
+    :param pandas.DataFrame lines: data frame of lines.
+    :param pandas.DataFrame substations: data frame of substations.
+    :return: (*list*) -- each entry is a 3-tuple:
+        index of the first connected component (int),
+        index of the second connected component (int),
+        a dictionary with keys containing ``start``, ``end`` and ``weight``, defines
+        the ``from substation ID``, ``to substation ID`` and the distance of the line.
+    """
+    cache_dir = os.path.join(os.path.dirname(__file__), "cache")
+    os.makedirs(cache_dir, exist_ok=True)
+    cache_key = repr(lines[["SUB_1_ID", "SUB_2_ID"]].to_numpy().tolist()) + repr(
+        substations[["LATITUDE", "LONGITUDE"]].to_numpy().tolist()
+    )
+    cache_hash = hashlib.md5(cache_key.encode("utf-8")).hexdigest()
+    try:
+        with open(os.path.join(cache_dir, f"mst_{cache_hash}.pkl"), "rb") as f:
+            print("Reading cached minimum spanning tree")
+            mst_edges = pickle.load(f)
+    except Exception:
+        print("No minimum spanning tree available, generating...")
+        _, mst_edges = connect_islands_with_minimum_cost(lines, substations)
+        with open(os.path.join(cache_dir, f"mst_{cache_hash}.pkl"), "wb") as f:
+            pickle.dump(mst_edges, f)
+    return mst_edges
+
+
 def build_transmission(method="sub2line", kwargs={"rounding": 3}):
     """Build transmission network
 
@@ -481,6 +599,27 @@ def build_transmission(method="sub2line", kwargs={"rounding": 3}):
             substations, hifld_lines, **kwargs
         )
 
+    # Connect all connected components
+    mst_edges = get_mst_edges(lines, substations)
+    new_lines = pd.DataFrame(
+        [{"SUB_1_ID": x[2]["start"], "SUB_2_ID": x[2]["end"]} for x in mst_edges]
+    )
+    new_lines = new_lines.assign(VOLTAGE=pd.NA, VOLT_CLASS="NOT AVAILABLE")
+    new_lines["COORDINATES"] = new_lines.apply(
+        lambda x: [
+            [
+                substations.loc[x.SUB_1_ID, "LATITUDE"],
+                substations.loc[x.SUB_1_ID, "LONGITUDE"],
+            ],
+            [
+                substations.loc[x.SUB_2_ID, "LATITUDE"],
+                substations.loc[x.SUB_2_ID, "LONGITUDE"],
+            ],
+        ],
+        axis=1,
+    )
+    lines = pd.concat([lines, new_lines])
+
     # Add voltages to lines with missing data
     augment_line_voltages(lines, substations)