diff --git a/docs/source/configtables/electricity.csv b/docs/source/configtables/electricity.csv index 7741f971..20b79ec6 100644 --- a/docs/source/configtables/electricity.csv +++ b/docs/source/configtables/electricity.csv @@ -52,6 +52,7 @@ demand_response:,,,Settings to activate and configure demand response ,,, imports:,,,Configure electric imports from regions outside of model scope -- enable,,``true`` or ``false``,Enable electric imports +-- representation,,``store`` or ``generator``,"How out-of-footprint supply is represented. ``store`` (default) puts a bottomless ``Store`` behind each external zone and PRICES the interface ``Link``. ``generator`` puts a priced import ``Generator`` (carrier ``unspecified_imports``, carrying ``co2_emissions``) on the external bus, leaves the interface link unpriced, and moves the CPUC contracted out-of-state units behind the boundary so their deliveries are metered against ``volume_limit`` and the interface capacity." -- costs,$/MWh,``wholesale`` or ```` or ``float``,Cost of electric imports from regions outside of model scope. ``Wholesale`` will use in monthly wholesales electric prices. ``float`` will assign a user specified value. ``carrier`` will take average marginal cost of the carrier. -- co2_emissions,CO2/MWh,``float``,CO2 emissions of electric imports from regions outside of model scope. -- capacity_limit,,``true`` or ``false``,Enable capacity limit for electric imports from regions outside of model scope diff --git a/workflow/repo_data/config/config.california.yaml b/workflow/repo_data/config/config.california.yaml index 4c0eb126..594e6a2b 100644 --- a/workflow/repo_data/config/config.california.yaml +++ b/workflow/repo_data/config/config.california.yaml @@ -55,6 +55,8 @@ electricity: profile: servm imports: enable: true + representation: generator # CPUC out-of-state contracts sit BEHIND the WECC boundary, so their deliveries are metered against the import cap + costs: 60 # flat $/MWh for unspecified imports; the 'wholesale' EIA series is retail-priced (see issue #807) volume_limit: 25 # % of annual CA load balancing_period: year exports: @@ -69,12 +71,6 @@ clustering: aggregation_strategies: generators: p_min_pu: 'capacity_weighted_average' # keep the UC feasibility invariant through clustering - # CPUC out-of-state units (REMOTE_PREFIX "R ") are attached to CA buses with - # p_nom_extendable=False, while in-state units of the same carrier are extendable - # because their carrier is in extendable_carriers.Generator. Clustering groups - # one-ports by (bus, carrier), so those two land in one group and PyPSA's default - # `consense` raises. Same situation `committable: any` already handles in the base. - p_nom_extendable: any # Hydro units disagree on up_time_before ({0, 1}). `max` keeps the integer # dtype (capacity_weighted_average would yield a float) and is the # UC-feasibility-safe direction: the cluster inherits the longest online diff --git a/workflow/repo_data/config/config.default.yaml b/workflow/repo_data/config/config.default.yaml index c23dd312..5214d8aa 100644 --- a/workflow/repo_data/config/config.default.yaml +++ b/workflow/repo_data/config/config.default.yaml @@ -228,6 +228,7 @@ electricity: # neighbors using historical EIA interchange data. imports: enable: false + representation: store # store = priced links off a bottomless Store; generator = priced generator (+ CPUC contracted units) behind an unpriced interface link costs: wholesale # wholesale | carrier | float — how imported energy is priced co2_emissions: 0.428 # tCO2/MWh assigned to imported energy capacity_limit: true # cap import power to historical maxima diff --git a/workflow/rules/build_electricity.smk b/workflow/rules/build_electricity.smk index cf542207..ebda8c08 100644 --- a/workflow/rules/build_electricity.smk +++ b/workflow/rules/build_electricity.smk @@ -804,6 +804,9 @@ rule add_electricity: planning_horizons=config["scenario"]["planning_horizons"], eia_api=config["api"]["eia"], remote_contracted=config["electricity"].get("remote_contracted_resources", {}), + imports_representation=config_provider( + "electricity", "imports", "representation", default="store" + ), input: unpack(dynamic_fuel_price_files), unpack(remote_contracted_resource_files), @@ -870,7 +873,11 @@ rule add_electricity: else [] ), output: - NETWORKS + "{interconnect}/elec_s{simpl}_l_pp.pkl", + network=NETWORKS + "{interconnect}/elec_s{simpl}_l_pp.pkl", + # CPUC contracted out-of-state units, serialized for attachment behind + # the model boundary when electricity.imports.representation == generator. + # Always declared; a sentinel is written otherwise. + remote_units=NETWORKS + "{interconnect}/remote_units_s{simpl}.pkl", log: LOGS + "{interconnect}/elec_s{simpl}_add_electricity.log", benchmark: @@ -1064,6 +1071,7 @@ rule add_extra_components: + "{interconnect}/regions_onshore_s{simpl}_{clusters}.geojson", flowgates=flowgates_for_extra_components, reeds_memberships="repo_data/ReEDS_Constraints/membership.csv", + remote_units=NETWORKS + "{interconnect}/remote_units_s{simpl}.pkl", co2_storage=( CO2 + "{interconnect}/co2_storage_s{simpl}_{clusters}.csv" if config["scenario"]["sector"] == "" and config["co2"]["storage"] is True @@ -1073,7 +1081,6 @@ rule add_extra_components: params: retirement=config["electricity"].get("retirement", "technical"), demand_response=config["electricity"].get("demand_response", {}), - trim_network=config_provider("model_topology", "trim", default=False), imports=config_provider("electricity", "imports", default={}), exports=config_provider("electricity", "exports", default={}), weather_year=config_provider("renewable_weather_years"), diff --git a/workflow/schemas/config.schema.yaml b/workflow/schemas/config.schema.yaml index faabe09a..0729b387 100644 --- a/workflow/schemas/config.schema.yaml +++ b/workflow/schemas/config.schema.yaml @@ -90,7 +90,6 @@ properties: transmission_network: {enum: [reeds, tamu]} topological_boundaries: {enum: [county, reeds_zone, state]} interface_transmission_limits: {type: boolean} - trim: {type: boolean} include: description: > Zone subset for footprint-scoped runs. null means "no filter"; see the @@ -285,7 +284,12 @@ properties: additionalProperties: false properties: enable: {type: boolean} - costs: {enum: [wholesale, carrier, float]} + representation: {enum: [store, generator]} + costs: + description: > + 'wholesale' (EIA price series), a carrier name (that fleet's + average marginal cost), or a flat number in $/MWh. + type: [string, number] co2_emissions: {type: number} capacity_limit: {type: boolean} volume_limit: {} @@ -296,7 +300,11 @@ properties: additionalProperties: false properties: enable: {type: boolean} - costs: {enum: [wholesale, carrier, float]} + costs: + description: > + 'wholesale' (EIA price series), a carrier name (that fleet's + average marginal cost), or a flat number in $/MWh. + type: [string, number] capacity_limit: {type: boolean} volume_limit: {} balancing_period: {enum: [day, week, month, year]} diff --git a/workflow/scripts/add_electricity.py b/workflow/scripts/add_electricity.py index c0c45b88..e57947f5 100755 --- a/workflow/scripts/add_electricity.py +++ b/workflow/scripts/add_electricity.py @@ -1301,6 +1301,40 @@ def attach_remote_contracted_resources( One row per ledger row with its disposition (``status``, ``carrier``, ``bus``, ``p_nom``), for logging and testing. """ + unit_df, summary = build_remote_contracted_units(n, plants_prefilter, remote_df, weights, costs, tech_map) + if unit_df.empty: + return summary + + dropped_vre = attach_remote_units(n, unit_df, costs, conventional_carriers, unit_commitment) + return _finalize_remote_summary(summary, dropped_vre) + + +def build_remote_contracted_units( + n: pypsa.Network, + plants_prefilter: pd.DataFrame, + remote_df: pd.DataFrame, + weights: pd.DataFrame, + costs: pd.DataFrame, + tech_map: pd.DataFrame | None = None, +) -> tuple[pd.DataFrame, pd.DataFrame]: + """Derive the CPUC contracted-unit table without touching the network. + + This is the half of :func:`attach_remote_contracted_resources` that only + needs data: it resolves each ledger row onto its live ``powerplants.csv`` + rows, derives capacity-weighted techno-economics, and picks the California + attachment bus (the region's max-LAF bus). In ``imports.representation: + generator`` the units are attached at an EXTERNAL bus instead, but the + California bus is still resolved here because remote VRE units borrow that + bus's capacity-factor profile. + + Returns + ------- + (unit_df, summary) + ``unit_df`` is indexed by ``REMOTE_PREFIX + cpuc_unit_name`` and carries + every attribute the ``n.add`` calls need, plus the physical ``state`` of + the constituent plants (used to place the unit behind the right external + zone). ``summary`` is the per-ledger-row disposition table. + """ region_bus = _servm_region_buses(weights, n) carrier_to_category, category_to_carrier = _carrier_category_maps(tech_map) @@ -1349,6 +1383,7 @@ def attach_remote_contracted_resources( bus = region_bus[row["servm_region"]] p_nom = float(min(record["capmax_mw"], constituents["p_nom"].sum())) attrs = _remote_unit_attributes(constituents, p_nom, carrier, costs) + state = _remote_unit_state(constituents) record.update( { @@ -1359,7 +1394,9 @@ def attach_remote_contracted_resources( }, ) records.append(record) - units.append({"name": REMOTE_PREFIX + name, "carrier": carrier, "bus": bus, "p_nom": p_nom, **attrs}) + units.append( + {"name": REMOTE_PREFIX + name, "carrier": carrier, "bus": bus, "p_nom": p_nom, "state": state, **attrs}, + ) logger.info( f"Remote contracted unit '{name}' -> bus '{bus}' ({row['servm_region']} max-LAF bus), carrier " @@ -1369,7 +1406,7 @@ def attach_remote_contracted_resources( summary = pd.DataFrame(records) if not units: logger.warning("Remote contracted resources enabled but no ledger row could be attached.") - return summary + return pd.DataFrame(), summary unit_df = pd.DataFrame(units).set_index("name") if unit_df.index.has_duplicates: @@ -1377,6 +1414,40 @@ def attach_remote_contracted_resources( f"Duplicate cpuc_unit_name(s) in the remote contracted-resource file: " f"{sorted(unit_df.index[unit_df.index.duplicated()])}", ) + return unit_df, summary + + +def _remote_unit_state(constituents: pd.DataFrame) -> str: + """Physical state of a contract's constituent plants (dominant by capacity). + + Hoover is the reason this is capacity-weighted rather than a single lookup: + the contract spans EIA 154 (NV) and 8902 (AZ), two halves of the same dam. + """ + if "state" not in constituents.columns: + return "" + states = constituents.dropna(subset=["state"]) + if states.empty: + return "" + return str(states.groupby("state")["p_nom"].sum().idxmax()) + + +def attach_remote_units( + n: pypsa.Network, + unit_df: pd.DataFrame, + costs: pd.DataFrame, + conventional_carriers: list, + unit_commitment: bool, + profiles: pd.DataFrame | None = None, +) -> list[str]: + """Add a derived contracted-unit table to the network at ``unit_df['bus']``. + + Shared by both import representations: in ``store`` mode the buses are + California buses and the profiles are derived here; in ``generator`` mode + ``external_regions`` rewrites ``bus`` to an external bus and supplies the + profiles that were borrowed at the pre-clustering stage. + + Returns the ledger names dropped for want of a VRE profile. + """ add_missing_carriers(n, sorted(set(unit_df["carrier"]))) batteries = unit_df[unit_df["carrier"] == "battery"] @@ -1384,9 +1455,48 @@ def attach_remote_contracted_resources( firm = unit_df.drop(index=batteries.index.union(vre.index)) _attach_remote_firm(n, firm, costs, conventional_carriers, unit_commitment) - dropped_vre = _attach_remote_vre(n, vre, costs) + dropped_vre = _attach_remote_vre(n, vre, costs, profiles=profiles) _attach_remote_batteries(n, batteries, costs) + return dropped_vre + +def build_remote_unit_bundle( + n: pypsa.Network, + unit_df: pd.DataFrame, + summary: pd.DataFrame, + costs: pd.DataFrame, + conventional_carriers: list, + unit_commitment: bool, +) -> dict: + """Serialize the contracted units instead of attaching them. + + Used by ``imports.representation: generator``, where the units belong at an + external bus that only exists later, in ``add_extra_components``. Everything + that stage cannot recompute is carried across: the derived unit table, the + CF profiles borrowed from the (pre-clustering) California buses, and the + cost table the ``n.add`` calls resolve ``capital_cost``/``lifetime`` from. + """ + vre = unit_df[unit_df["carrier"].isin(VRE_PROFILE_CARRIERS)] + profiles, dropped_vre = collect_remote_vre_profiles(n, vre) + unit_df = unit_df.drop(index=[REMOTE_PREFIX + name for name in dropped_vre], errors="ignore") + + bundle = { + "units": unit_df, + "vre_profiles": pd.DataFrame(profiles) if profiles else pd.DataFrame(index=n.snapshots), + "costs": costs, + "conventional_carriers": list(conventional_carriers), + "unit_commitment": bool(unit_commitment), + "summary": _finalize_remote_summary(summary, dropped_vre), + } + logger.info( + f"Serialized {len(unit_df)} remote contracted unit(s) ({unit_df['p_nom'].sum():.1f} MW) for attachment " + "behind the model boundary (imports.representation: generator).", + ) + return bundle + + +def _finalize_remote_summary(summary: pd.DataFrame, dropped_vre: list[str]) -> pd.DataFrame: + """Fold the dropped-VRE outcome into the summary and log the totals.""" if dropped_vre: summary.loc[summary["cpuc_unit_name"].isin(dropped_vre), "status"] = "skipped_no_profile" summary.loc[summary["cpuc_unit_name"].isin(dropped_vre), "p_nom"] = 0.0 @@ -1470,16 +1580,14 @@ def _attach_remote_firm( _apply_remote_seasonal_derates(n, firm[["summer_derate", "winter_derate"]]) -def _attach_remote_vre(n: pypsa.Network, vre: pd.DataFrame, costs: pd.DataFrame) -> list[str]: - """Attach remote wind/solar contracts, borrowing the attachment bus's CF profile. +def collect_remote_vre_profiles(n: pypsa.Network, vre: pd.DataFrame) -> tuple[dict[str, pd.Series], list[str]]: + """Borrow a CF profile per remote VRE unit from its attachment bus. - Returns the ledger names that had to be dropped for want of any profile. + Returns ``(profiles, dropped)``; ``dropped`` holds the ledger names for which + the network carries no profile of that carrier at all. """ + profiles: dict[str, pd.Series] = {} dropped: list[str] = [] - if vre.empty: - return dropped - - profiles = {} for name, unit in vre.iterrows(): profile = _remote_vre_profile(n, unit["bus"], unit["carrier"]) if profile is None: @@ -1490,12 +1598,39 @@ def _attach_remote_vre(n: pypsa.Network, vre: pd.DataFrame, costs: pd.DataFrame) dropped.append(name[len(REMOTE_PREFIX) :]) continue profiles[name] = profile + return profiles, dropped + + +def _attach_remote_vre( + n: pypsa.Network, + vre: pd.DataFrame, + costs: pd.DataFrame, + profiles: pd.DataFrame | None = None, +) -> list[str]: + """Attach remote wind/solar contracts, borrowing the attachment bus's CF profile. + + ``profiles`` short-circuits the borrowing: in ``imports.representation: + generator`` the units land on an external bus that has no profile of its + own, so the profiles borrowed before clustering are carried in on the + serialized bundle instead. + + Returns the ledger names that had to be dropped for want of any profile. + """ + dropped: list[str] = [] + if vre.empty: + return dropped + + if profiles is None: + profiles, dropped = collect_remote_vre_profiles(n, vre) + else: + profiles = {name: profiles[name] for name in vre.index if name in profiles.columns} vre = vre.loc[list(profiles)] if vre.empty: return dropped p_max_pu = pd.DataFrame(profiles).reindex(columns=vre.index) + p_max_pu = p_max_pu.set_axis(n.snapshots) n.add( "Generator", vre.index, @@ -1825,6 +1960,7 @@ def main(snakemake): # (small) slice of the fleet the CPUC ledger points at BEFORE that filter runs; # attach_remote_contracted_resources adds them back at CA buses later. remote_contracted = dict(getattr(params, "remote_contracted", None) or {}) + imports_representation = getattr(params, "imports_representation", "store") or "store" remote_df = None plants_prefilter = None if remote_contracted.get("enable", False): @@ -1924,18 +2060,42 @@ def main(snakemake): # Runs last among the attach_* steps: remote VRE contracts copy the # capacity-factor profile of their California attachment bus, which only # exists once attach_wind_and_solar has run. + # + # `imports.representation: generator` moves these units BEHIND the model + # boundary, so they must not be added here — add_extra_components attaches + # them at the external import buses instead. The bundle output is declared + # unconditionally by the rule, so a sentinel is written when there is + # nothing to hand over. + remote_bundle = None if remote_df is not None: - attach_remote_contracted_resources( + unit_df, summary = build_remote_contracted_units( n, plants_prefilter, remote_df, pd.read_csv(snakemake.input["servm_load_weights"]), costs, - conventional_carriers, - extendable_carriers, tech_map=pd.read_csv(snakemake.input["servm_tech_map"]), - unit_commitment=params.conventional["unit_commitment"], ) + if unit_df.empty: + pass + elif imports_representation == "generator": + remote_bundle = build_remote_unit_bundle( + n, + unit_df, + summary, + costs, + conventional_carriers, + params.conventional["unit_commitment"], + ) + else: + dropped_vre = attach_remote_units( + n, + unit_df, + costs, + conventional_carriers, + params.conventional["unit_commitment"], + ) + _finalize_remote_summary(summary, dropped_vre) update_p_nom_max(n) @@ -2010,7 +2170,7 @@ def main(snakemake): axis=1, ) - output_folder = os.path.dirname(snakemake.output[0]) + "/base_network" + output_folder = os.path.dirname(snakemake.output.network) + "/base_network" export_network_for_gis_mapping(n, output_folder) clean_bus_data(n) @@ -2018,8 +2178,12 @@ def main(snakemake): n.meta = snakemake.config log_network_schema(n, stage="exit", baseline=schema_entry) - # n.export_to_netcdf(snakemake.output[0]) - pickle.dump(n, open(snakemake.output[0], "wb")) + # n.export_to_netcdf(snakemake.output.network) + pickle.dump(n, open(snakemake.output.network, "wb")) + + # Always written, even when it holds nothing: snakemake requires every + # declared output to exist. + pickle.dump(remote_bundle, open(snakemake.output.remote_units, "wb")) if __name__ == "__main__": diff --git a/workflow/scripts/add_extra_components.py b/workflow/scripts/add_extra_components.py index 061ddc93..b167d5f2 100644 --- a/workflow/scripts/add_extra_components.py +++ b/workflow/scripts/add_extra_components.py @@ -9,8 +9,13 @@ from _helpers import calculate_annuity, configure_logging, load_costs, log_network_schema from add_electricity import add_missing_carriers from constants import HOURS_PER_YEAR -from eia import FuelCosts -from opts._helpers import get_region_buses +from external_regions import ( + add_external_regions, + convert_flowgates_to_state, + format_flowgates_for_imports_exports, + load_remote_unit_bundle, + resolve_trade_costs, +) from shapely.geometry import Point idx = pd.IndexSlice @@ -696,370 +701,6 @@ def add_demand_response( ) -def trim_network(n, trim_topology): - """ - Trim_network splits the network into two parts: - - The internal network, which is the network within the specified zones. - - The external network, which is the network outside the specified zones. - - The internal network is retained and unchanged. While the external network components are removed. The external buses which are directly connected to the internal network are aggregated to the `nerc_reg` value of their buses. - The only generators kept are the OCGTs at the external buses, which are set to non-extendable. - - The external OCGT generators are set to the carrier name `imports` and retain the same emissions intensity. - - """ - retain_zones = trim_topology["zone"] - internal_buses = get_region_buses(n, retain_zones) - if internal_buses.empty: - logger.warning("No internal buses found, skipping trim_network") - return None - - # Get all lines and links connected to internal buses - retain_lines = n.lines[n.lines.bus0.isin(internal_buses.index) | n.lines.bus1.isin(internal_buses.index)] - retain_links = n.links[n.links.bus0.isin(internal_buses.index) | n.links.bus1.isin(internal_buses.index)] - - # Find buses to remove (those not connected to internal network) - buses_to_remove = n.buses[ - ~n.buses.index.isin(retain_lines.bus0) - & ~n.buses.index.isin(retain_lines.bus1) - & ~n.buses.index.isin(retain_links.bus0) - & ~n.buses.index.isin(retain_links.bus1) - ] - - # Find external buses to keep (connected to internal network but not internal) - external_buses_to_keep = n.buses.loc[ - ~n.buses.index.isin(buses_to_remove.index) & ~n.buses.index.isin(internal_buses.index) - ] - - # Remove components at buses that are being removed - for c in n.one_port_components: - component = n.components[c].static - rm = component[component.bus.isin(buses_to_remove.index)] - if not rm.empty: - n.remove(c, rm.index) - - # Remove lines and links at buses being removed - for c in ["Line", "Link"]: - component = n.components[c].static - rm = component[~component.bus0.isin(internal_buses.index) & ~component.bus1.isin(internal_buses.index)] - if not rm.empty: - n.remove(c, rm.index) - - # Remove the buses - n.remove("Bus", buses_to_remove.index) - - # Get OCGT generators and calculate average marginal cost - ocgt_gens = n.generators[n.generators.carrier == "OCGT"] - avg_marginal_cost = n.get_switchable_as_dense("Generator", "marginal_cost").loc[:, ocgt_gens.index].mean().mean() - n.add("Carrier", "imports", co2_emissions=0.428, nice_name="imports") - - # remove existing oneport components at bus - for c in n.one_port_components: - component = n.components[c].static - rm = component[component.bus.isin(external_buses_to_keep.index)] - if not rm.empty: - logger.info(f"Removing {c} at external buses {external_buses_to_keep.index} with components {rm.index}") - n.remove(c, rm.index) - - # Handle external buses and their generators - for bus in external_buses_to_keep.index: - # Create new import generator - bus_name = n.buses.loc[bus].name - n.add( - "Generator", - f"import_{bus_name}", - bus=bus, - carrier="imports", - p_nom=1e4, - p_nom_extendable=False, - marginal_cost=avg_marginal_cost, - efficiency=1, - build_year=n.investment_periods[0], - lifetime=100, - ) - - # Change location names of external buses, append imports to the ['reeds_state', 'reeds_zone', 'reeds_ba', 'interconnect', 'trans_reg', 'trans_grp'] - n.buses.loc[bus, "reeds_state"] = f"imports_{n.buses.loc[bus, 'reeds_state']}" - n.buses.loc[bus, "reeds_zone"] = f"imports_{n.buses.loc[bus, 'reeds_zone']}" - n.buses.loc[bus, "reeds_ba"] = f"imports_{n.buses.loc[bus, 'reeds_ba']}" - n.buses.loc[bus, "interconnect"] = f"imports_{n.buses.loc[bus, 'interconnect']}" - n.buses.loc[bus, "trans_reg"] = f"imports_{n.buses.loc[bus, 'trans_reg']}" - n.buses.loc[bus, "trans_grp"] = f"imports_{n.buses.loc[bus, 'trans_grp']}" - - # Set all links and lines connected to the bus as non-extendable - for c in ["Line", "Link"]: - attr_name = "p_nom_extendable" if c == "Link" else "s_nom_extendable" - component = n.components[c].static - mask = (component.bus0 == bus) | (component.bus1 == bus) - if mask.any(): - component.loc[mask, attr_name] = False - n.components[c].static.update(component) - - # Remove the links which have "exp" in the name and are connected to the external buses - links_to_remove = n.links[ - n.links.index.str.contains("exp") - & (n.links.bus0.isin(external_buses_to_keep.index) | n.links.bus1.isin(external_buses_to_keep.index)) - ] - n.remove("Link", links_to_remove.index) - - # Update network topology - n.determine_network_topology() - - -def calc_import_export_costs(n: pypsa.Network, carrier: str) -> float: - """Calculates the average marginal cost for a given carrier.""" - gens = n.generators[n.generators.carrier == carrier] - component = "Generator" - if gens.empty: - gens = n.links[n.links.carrier == carrier] - component = "Link" - if gens.empty: - raise ValueError(f"No generators or links found for carrier to calculate imports/exports costs: {carrier}") - costs = n.get_switchable_as_dense(component, "marginal_cost").loc[:, gens.index].mean().mean() - if costs <= 0.01: - raise ValueError( - f"Average marginal cost for {carrier} is less than or equal to 0.01. Check the fuel costs configuration.", - ) - return costs - - -def load_import_export_costs(eia_api: str, year: int) -> pd.DataFrame: - """Loads fuel costs from EIA.""" - return FuelCosts(fuel="electricity", year=year, api=eia_api).get_data() - - -def format_import_export_costs(n: pypsa.Network, fuel_costs: pd.DataFrame) -> pd.DataFrame: - """Formats fuel costs for BA mappings.""" - df = fuel_costs.copy() - data = [] - - buses = n.buses.copy() - - region_mapping = buses.set_index("country")["reeds_state"].to_dict() - for region, state in region_mapping.items(): - for period in df.index.unique(): - temp = df[(df.index == period) & (df.state == state)] - value = temp.value.mean() - data.append([period, region, value, "usd/mwh"]) - formatted = pd.DataFrame(data, columns=["period", "zone", "value", "units"]).set_index("period") - return formatted[~formatted.value.isna()] # regions outside of model scope - - -def format_flowgates_for_imports_exports(n: pypsa.Network, flowgates: pd.DataFrame, zone_col: str) -> pd.DataFrame: - """Formats flowgates for zone mappings.""" - zones_in_model = n.buses[zone_col].unique() - df = flowgates.copy() - - # only keep flowgates that connect inside to outside model scope - df = df[df.r.isin(zones_in_model) ^ df.rr.isin(zones_in_model)] - - # reformat to sinlge value column for easier addition to network - data = [] - for _, row in df.iterrows(): - if row.MW_f0 > 0: - data.append([row.r, row.rr, row.MW_f0]) - if row.MW_r0 > 0: - data.append([row.rr, row.r, row.MW_r0]) - - return pd.DataFrame(data, columns=["r", "rr", "value"]) - - -def convert_flowgates_to_state(flowgates: pd.DataFrame, membership: pd.DataFrame) -> pd.DataFrame: - """Converts flowgates to state level.""" - mbshp = membership.set_index("ba") - df = flowgates.copy() - - df["s"] = df.r.map(mbshp["st"]) - df["ss"] = df.rr.map(mbshp["st"]) - df = df.drop(columns=["r", "rr"]) - df = df.rename(columns={"s": "r", "ss": "rr"}) - return df - - -def add_elec_imports_exports( - n: pypsa.Network, - direction: str, - flowgates: pd.DataFrame, - fuel_costs: pd.DataFrame | float, - co2_emissions: float = 0, - zone_col: str = "reeds_zone", -): - """Add electricity imports and exports to the network. - - These are capacity constrianed links to/from states outside the model spatial scope. - """ - - def _get_regions_2_add(n: pypsa.Network, flowgates: pd.DataFrame, zone_col: str) -> list[str]: - """Gets regions to add import and export buses to.""" - unique_regions = set(flowgates.r.unique()) | set(flowgates.rr.unique()) - return [x for x in unique_regions if x not in n.buses[zone_col].unique()] - - def _add_import_export_carriers(n: pypsa.Network, direction: str, co2_emissions: float | None = None) -> None: - """Adds import and export carriers to the network.""" - if direction == "imports": - co2_emissions = 0 if not co2_emissions else co2_emissions - n.add("Carrier", "imports", co2_emissions=co2_emissions, nice_name="Imports") - elif direction == "exports": - n.add("Carrier", "exports", co2_emissions=0, nice_name="Exports") - else: - raise ValueError(f"direction must be either imports or exports; received: {direction}") - - def _add_import_export_buses(n: pypsa.Network, regions_2_add: list[str], direction: str) -> None: - """Adds import and export buses to the network.""" - if direction == "imports": - suffix = "_imports" - carrier = "imports" - elif direction == "exports": - suffix = "_exports" - carrier = "exports" - else: - raise ValueError(f"direction must be either imports or exports; received: {direction}") - - # cant add in the reeds_state, reeds_zone, reeds_ba, interconnect, trans_reg, trans_grp - # because this information has already been filtered out of the network - - n.add( - "Bus", - regions_2_add, - suffix=suffix, - carrier=carrier, - country=regions_2_add, - ) - - def _add_import_export_stores(n: pypsa.Network, regions_2_add: list[str], direction: str) -> None: - """Adds import and export stores to the network.""" - if direction == "imports": - n.add( - "Store", - regions_2_add, - bus=[f"{x}_imports" for x in regions_2_add], - suffix="_imports", - carrier="imports", - e_nom=0, - e_nom_extendable=True, - capital_cost=0, - e_nom_min=0, - e_nom_max=1e9, - e_min_pu=-1, - e_max_pu=0, - e_cyclic_per_period=False, - marginal_cost=0, - ) - elif direction == "exports": - n.add( - "Store", - regions_2_add, - bus=[f"{x}_exports" for x in regions_2_add], - suffix="_exports", - carrier="exports", - e_nom_extendable=True, - marginal_cost=0, - e_nom=0, - e_nom_max=1e9, - e_min=0, - e_min_pu=0, - e_max_pu=1, - ) - else: - raise ValueError(f"direction must be either imports or exports; received: {direction}") - - def _build_cost_timeseries(n: pypsa.Network, costs: pd.DataFrame, zone: str) -> pd.Series: - """Builds a cost timeseries for a given state.""" - timesteps = n.snapshots.get_level_values("timestep") - years = n.investment_periods - cost_by_zone = costs[costs.zone == zone].drop(columns=["zone", "units"]) - dfs = [] - for year in years: - df = cost_by_zone.copy() - df.index = pd.to_datetime(df.index).map(lambda x: x.replace(year=year)) - df = df.resample("h").ffill().reindex(timesteps).ffill() - df["year"] = year - df = df.set_index(["year", df.index]) # df.index is timestep - dfs.append(df) - df = pd.concat(dfs) - return df.reindex(n.snapshots) - - def _add_import_export_links( - n: pypsa.Network, - flowgates: pd.DataFrame, - fuel_costs: pd.DataFrame | float | str, - direction: str, - zone_col: str = "reeds_zone", - ) -> None: - """Adds import and export links to the network.""" - costs = {} - zones_in_model = n.buses[zone_col].dropna().unique() - - for _, row in flowgates.iterrows(): - zone_inside = row.r if row.r in zones_in_model else row.rr - zone_outside = row.r if row.r not in zones_in_model else row.rr - - # extremely crude caching for generating cost timeseries :| - # keyed by the INSIDE zone — checking the outside zone here skipped - # the write whenever an earlier row's inside zone happened to match, - # leaving costs[zone_inside] unset (KeyError on county networks). - if zone_inside not in costs: - if isinstance(fuel_costs, float | int): - costs[zone_inside] = fuel_costs - elif isinstance(fuel_costs, pd.DataFrame): - costs[zone_inside] = _build_cost_timeseries(n, fuel_costs, zone_inside) - else: - costs[zone_inside] = 0 - - marginal_cost = costs[zone_inside] - - capacity = row.value - - """Structre of flowgates is given by: - - r rr value - 0 p6 p8 488.117 - 1 p8 p6 378.458 - 2 p6 p9 4800.000 - ... - """ - - if direction == "imports": - if row.r == zone_inside: # originating at r is exports (ie r -> rr) - continue - name = f"{zone_inside}_{zone_outside}_imports" - bus0 = f"{zone_outside}_imports" - bus1 = zone_inside - carrier = "imports" - else: - if row.r == zone_outside: # originating at rr is exports (ie rr -> r) - continue - name = f"{zone_inside}_{zone_outside}_exports" - bus0 = zone_inside - bus1 = f"{zone_outside}_exports" - carrier = "exports" - if isinstance(marginal_cost, pd.Series): - marginal_cost = marginal_cost.mul(-1) # constraint will limit exports - - mc = marginal_cost.value if isinstance(marginal_cost, pd.DataFrame) else marginal_cost - - n.add( - "Link", - name, - bus0=bus0, - bus1=bus1, - carrier=carrier, - p_nom_extendable=False, - p_min_pu=0, - p_max_pu=1, - marginal_cost=mc, - p_nom=capacity, - ) - - assert direction in ["imports", "exports"], f"direction must be either imports or exports; received: {direction}" - - regions_2_add = _get_regions_2_add(n, flowgates, zone_col) - _add_import_export_carriers(n, direction, co2_emissions) - _add_import_export_buses(n, regions_2_add, direction) - _add_import_export_stores(n, regions_2_add, direction) - _add_import_export_links(n, flowgates, fuel_costs, direction, zone_col) - - def add_co2_storage(n: pypsa.Network, config: dict, co2_storage_csv: str, costs: pd.DataFrame, sector: bool): """Adds node level CO2 (underground) storage.""" # get node level CO2 (underground) storage potential and cost from CSV file @@ -1556,23 +1197,16 @@ def main(snakemake) -> None: if dr_config: add_demand_response(n, dr_config) - trim_network_config = snakemake.params.trim_network imports_config = snakemake.params.imports exports_config = snakemake.params.exports - - assert not ( - snakemake.params.trim_network and (imports_config.get("enable", False) or exports_config.get("enable", False)) - ), "trim_network and imports/exports cannot be used together" - - if snakemake.params.trim_network: - trim_network(n, trim_network_config) + representation = imports_config.get("representation", "store") if snakemake.params.transmission_network == "reeds": # flowgates to limit the capacity (removed later if configured capacity limit is inf) flowgates = pd.read_csv(snakemake.input.flowgates) + membership = pd.read_csv(snakemake.input.reeds_memberships) if snakemake.params.topological_boundaries == "state": zone_col = "reeds_state" - membership = pd.read_csv(snakemake.input.reeds_memberships) flowgates = convert_flowgates_to_state(flowgates, membership) flowgates = format_flowgates_for_imports_exports(n, flowgates, zone_col) flowgates = flowgates.groupby(["r", "rr"], as_index=False).sum() @@ -1597,22 +1231,24 @@ def main(snakemake) -> None: if not imports_config.get("capacity_limit", True): import_flowgates["value"] = np.inf - import_costs = imports_config.get("costs", False) + fuel_costs = resolve_trade_costs(n, imports_config, "imports", snakemake.params.eia_api, year) - if isinstance(import_costs, float | int): # user defined value - fuel_costs = import_costs - elif isinstance(import_costs, str): # 'wholesale' or name of carrier - if import_costs == "wholesale": - fuel_costs = load_import_export_costs(snakemake.params.eia_api, year) - fuel_costs = format_import_export_costs(n, fuel_costs) - else: - fuel_costs = calc_import_export_costs(n, import_costs) - else: - raise ValueError( - f"'imports.costs' must be 'wholesale', name of a carrier, or a float/int. Received: {import_costs}", - ) + # Only `generator` mode places the CPUC contracted units behind the + # boundary; in `store` mode add_electricity has already attached them + # and the bundle is a sentinel. + remote_bundle = load_remote_unit_bundle(snakemake.input.remote_units) if representation == "generator" else None - add_elec_imports_exports(n, "imports", import_flowgates, fuel_costs, co2_emissions, zone_col) + add_external_regions( + n, + "imports", + representation, + import_flowgates, + fuel_costs, + co2_emissions, + zone_col, + remote_bundle=remote_bundle, + membership=membership, + ) # Electricity exports configuration if exports_config.get("enable", False) and snakemake.params.transmission_network == "reeds": @@ -1627,25 +1263,17 @@ def main(snakemake) -> None: if not exports_config.get("capacity_limit", True): export_flowgates["value"] = np.inf - export_costs = exports_config.get("costs", False) + fuel_costs = resolve_trade_costs(n, exports_config, "exports", snakemake.params.eia_api, year) - if isinstance(export_costs, float | int): # user defined value - fuel_costs = export_costs - fuel_costs *= -1 # make money by exporting - elif isinstance(export_costs, str): # 'wholesale' or name of carrier - if export_costs == "wholesale": - fuel_costs = load_import_export_costs(snakemake.params.eia_api, year) - fuel_costs = format_import_export_costs(n, fuel_costs) - fuel_costs["value"] = fuel_costs.value.mul(-1) # make money by exporting - else: - fuel_costs = calc_import_export_costs(n, export_costs) - fuel_costs *= -1 # make money by exporting - else: - raise ValueError( - f"'exports.costs' must be 'wholesale', name of a carrier, or a float/int. Received: {export_costs}", - ) - - add_elec_imports_exports(n, "exports", export_flowgates, fuel_costs, co2_emissions, zone_col) + add_external_regions( + n, + "exports", + representation, + export_flowgates, + fuel_costs, + co2_emissions, + zone_col, + ) if snakemake.config["scenario"]["sector"] == "E": co2_storage = snakemake.config.get("co2", {}).get("storage", False) diff --git a/workflow/scripts/build_fuel_prices.py b/workflow/scripts/build_fuel_prices.py index dd29aa58..482ffe09 100644 --- a/workflow/scripts/build_fuel_prices.py +++ b/workflow/scripts/build_fuel_prices.py @@ -61,24 +61,43 @@ def make_hourly(df: pd.DataFrame) -> pd.DataFrame: ### -def get_state_ng_power_prices(sns: pd.date_range, eia_api: str) -> pd.DataFrame: - df = ( - eia.FuelCosts("gas", sns.year[0], eia_api, industry="power").get_data( - pivot=True, +def _clamped_fuel_costs( + fuel: str, + sns: pd.date_range, + eia_api: str, + floor: int, +) -> pd.DataFrame: + """ + Fetch EIA fuel costs, clamping the query year to the earliest with data. + + EIA electric-power gas prices begin 2002-01 and coal shipment receipts + begin 2008; earlier years return empty payloads that crash format_data. + For earlier snapshot years fetch the floor year and shift the month-start + index back onto the snapshot year (a shift, not a replace, because the + inclusive API end bound leaves a spillover January of the next year). + """ + year = sns.year[0] + data_year = max(year, floor) + df = eia.FuelCosts(fuel, data_year, eia_api, industry="power").get_data( + pivot=True, + ) + if data_year != year: + logger.warning( + f"No EIA {fuel} power prices before {floor}; using {data_year} prices for snapshot year {year}", + ) + df.index = df.index.map( + lambda ts: ts.replace(year=ts.year - (data_year - year)), ) - * 1000 - / const.NG_MWH_2_MMCF - ) # $/MCF -> $/MWh + return df + + +def get_state_ng_power_prices(sns: pd.date_range, eia_api: str) -> pd.DataFrame: + df = _clamped_fuel_costs("gas", sns, eia_api, floor=2002) * 1000 / const.NG_MWH_2_MMCF # $/MCF -> $/MWh return make_hourly(df) def get_state_coal_power_prices(sns: pd.date_range, eia_api: str) -> pd.DataFrame: - eia_coal = ( - eia.FuelCosts("coal", sns.year[0], eia_api, industry="power").get_data( - pivot=True, - ) - * const.COAL_dol_ton_2_MWHthermal - ) + eia_coal = _clamped_fuel_costs("coal", sns, eia_api, floor=2008) * const.COAL_dol_ton_2_MWHthermal return make_hourly(eia_coal) diff --git a/workflow/scripts/external_regions.py b/workflow/scripts/external_regions.py new file mode 100644 index 00000000..6038258d --- /dev/null +++ b/workflow/scripts/external_regions.py @@ -0,0 +1,767 @@ +"""External (out-of-footprint) regions: imports, exports and contracted units. + +A regionally scoped run (California, say) cuts the synchronous grid at a +political boundary. Everything behind that cut still serves load inside the +footprint, and this module is the single place where it is represented. Two +representations are available, selected by ``electricity.imports.representation``: + +``store`` (default, unchanged behaviour) + Per external flowgate zone, a ``{zone}_imports`` bus carrying a bottomless + ``Store`` (``e_nom_max`` 1e9, ``e_min_pu`` -1) and one-way ``Link`` s into + the internal zone buses. Each link is rated at the NARIS flowgate capacity + and PRICED at the import price (a ``wholesale`` EIA timeseries, a carrier's + average marginal cost, or a flat float). Emissions are carried by the + ``imports`` carrier that the Store belongs to. Exports mirror this with an + absorbing Store and negatively-priced links. + +``generator`` + The external zone is modelled as a place with generation rather than as an + infinite energy tank: + + * the ``{zone}_imports`` bus carries a generic import ``Generator`` + (carrier ``unspecified_imports``) sized at the zone's total inbound + interface capacity and priced with the same machinery as ``store`` mode; + * California's CPUC-contracted out-of-state units (Palo Verde, Intermountain, + Hoover, Apex, the AZ/NV solar and battery contracts) are attached at that + same external bus instead of at a California bus, i.e. BEHIND the boundary; + * the ``Link`` s into the footprint are UNPRICED and rated at the flowgate + capacity, so the price sits on the generator and the link is pure transfer + capacity. + + The point of the second mode is accounting: because deliveries now traverse + a carrier-``imports`` link, they are counted by + :func:`opts.interchange.add_interchange_constraints` against the import + volume cap, and bounded by the interface capacity, exactly like generic + imports. In ``store`` mode the contracted units sit inside the footprint and + look like in-state generation, bypassing both limits. + +Both modes keep the link carriers ``imports`` / ``exports`` untouched, which is +what ``opts/interchange.py`` and ``opts/interfaces.py`` key off. + +External-bus naming +------------------- +Both modes use SEPARATE ``{zone}_imports`` and ``{zone}_exports`` buses rather +than one shared ``{zone}_external`` bus. Beyond keeping ``store`` mode +byte-identical, the separation is load-bearing in ``generator`` mode: a shared +bus would let the priced import generator (and the contracted units) sell +straight into the export sink, collecting the export price for free whenever the +export price exceeds the import cost — a pure arbitrage loop with no physical +meaning. + +Export pricing +-------------- +Export pricing is IDENTICAL in both modes: the negative price stays on the +export ``Link`` and the absorbing ``Store`` behind it is free. On the export +side there is no double-charge to worry about, because the only way energy can +reach the export Store is through exactly one export link — nothing else injects +into a ``{zone}_exports`` bus — so the negative price is earned exactly once per +exported MWh. That is what makes it safe to leave the export half of the +construction untouched by the representation switch, and it is only safe because +the import generator lives on a different bus (see above). + +CO2 +--- +In ``store`` mode the ``imports`` carrier carries ``imports.co2_emissions`` and +the Store's withdrawal is what the global CO2 constraint sees. In ``generator`` +mode there is no import Store, so the emission factor moves onto the +``unspecified_imports`` carrier of the generic import generator and ``imports`` +is set to zero. PyPSA attributes primary-energy emissions to generators and +stores, never to links, so the carrier-``imports`` links never double-count. +""" + +import logging +import re + +import dill +import pandas as pd +import pypsa +from add_electricity import add_missing_carriers, attach_remote_units +from eia import FuelCosts + +logger = logging.getLogger(__name__) + +REPRESENTATIONS = ("store", "generator") + +#: Carrier of the generic import generator in ``generator`` mode. Deliberately +#: NOT ``imports``: that carrier is reserved for the transfer links which +#: ``opts/interchange.py`` sums over. +GENERIC_IMPORT_CARRIER = "unspecified_imports" + + +# --------------------------------------------------------------------------- +# Flowgate formatting +# --------------------------------------------------------------------------- + + +def format_flowgates_for_imports_exports(n: pypsa.Network, flowgates: pd.DataFrame, zone_col: str) -> pd.DataFrame: + """Formats flowgates for zone mappings.""" + zones_in_model = n.buses[zone_col].unique() + df = flowgates.copy() + + # only keep flowgates that connect inside to outside model scope + df = df[df.r.isin(zones_in_model) ^ df.rr.isin(zones_in_model)] + + # reformat to sinlge value column for easier addition to network + data = [] + for _, row in df.iterrows(): + if row.MW_f0 > 0: + data.append([row.r, row.rr, row.MW_f0]) + if row.MW_r0 > 0: + data.append([row.rr, row.r, row.MW_r0]) + + return pd.DataFrame(data, columns=["r", "rr", "value"]) + + +def convert_flowgates_to_state(flowgates: pd.DataFrame, membership: pd.DataFrame) -> pd.DataFrame: + """Converts flowgates to state level.""" + mbshp = membership.set_index("ba") + df = flowgates.copy() + + df["s"] = df.r.map(mbshp["st"]) + df["ss"] = df.rr.map(mbshp["st"]) + df = df.drop(columns=["r", "rr"]) + df = df.rename(columns={"s": "r", "ss": "rr"}) + return df + + +# --------------------------------------------------------------------------- +# Import / export pricing +# --------------------------------------------------------------------------- + + +def calc_import_export_costs(n: pypsa.Network, carrier: str) -> float: + """Calculates the average marginal cost for a given carrier.""" + gens = n.generators[n.generators.carrier == carrier] + component = "Generator" + if gens.empty: + gens = n.links[n.links.carrier == carrier] + component = "Link" + if gens.empty: + raise ValueError(f"No generators or links found for carrier to calculate imports/exports costs: {carrier}") + costs = n.get_switchable_as_dense(component, "marginal_cost").loc[:, gens.index].mean().mean() + if costs <= 0.01: + raise ValueError( + f"Average marginal cost for {carrier} is less than or equal to 0.01. Check the fuel costs configuration.", + ) + return costs + + +def load_import_export_costs(eia_api: str, year: int) -> pd.DataFrame: + """Loads fuel costs from EIA.""" + # EIA retail-sales electricity prices begin 2001; earlier years return an + # empty payload that crashes format_data. No date shifting is needed: the + # downstream _build_cost_timeseries relabels the index onto the network's + # investment periods anyway. + data_year = max(year, 2001) + if data_year != year: + logger.warning( + f"No EIA electricity prices before 2001; using {data_year} prices for year {year}", + ) + return FuelCosts(fuel="electricity", year=data_year, api=eia_api).get_data() + + +def format_import_export_costs(n: pypsa.Network, fuel_costs: pd.DataFrame) -> pd.DataFrame: + """Formats fuel costs for BA mappings.""" + df = fuel_costs.copy() + data = [] + + buses = n.buses.copy() + + region_mapping = buses.set_index("country")["reeds_state"].to_dict() + for region, state in region_mapping.items(): + for period in df.index.unique(): + temp = df[(df.index == period) & (df.state == state)] + value = temp.value.mean() + data.append([period, region, value, "usd/mwh"]) + formatted = pd.DataFrame(data, columns=["period", "zone", "value", "units"]).set_index("period") + return formatted[~formatted.value.isna()] # regions outside of model scope + + +def resolve_trade_costs( + n: pypsa.Network, + trade_config: dict, + direction: str, + eia_api: str, + year: int, +) -> pd.DataFrame | float: + """Resolve ``imports.costs`` / ``exports.costs`` into a price the network can use. + + ``wholesale`` pulls the monthly EIA electricity price per state, a carrier + name averages that carrier's marginal cost, and a float is taken at face + value. Export prices are negated: exporting earns money. + """ + if direction not in ("imports", "exports"): + raise ValueError(f"direction must be either imports or exports; received: {direction}") + + sign = 1 if direction == "imports" else -1 + costs = trade_config.get("costs", False) + + if isinstance(costs, float | int): # user defined value + return costs * sign + if isinstance(costs, str): # 'wholesale' or name of carrier + if costs == "wholesale": + fuel_costs = load_import_export_costs(eia_api, year) + fuel_costs = format_import_export_costs(n, fuel_costs) + if sign < 0: + fuel_costs["value"] = fuel_costs.value.mul(sign) # make money by exporting + return fuel_costs + return calc_import_export_costs(n, costs) * sign + raise ValueError( + f"'{direction}.costs' must be 'wholesale', name of a carrier, or a float/int. Received: {costs}", + ) + + +def _build_cost_timeseries(n: pypsa.Network, costs: pd.DataFrame, zone: str) -> pd.Series: + """Builds a cost timeseries for a given state.""" + timesteps = n.snapshots.get_level_values("timestep") + years = n.investment_periods + cost_by_zone = costs[costs.zone == zone].drop(columns=["zone", "units"]) + dfs = [] + for year in years: + df = cost_by_zone.copy() + df.index = pd.to_datetime(df.index).map(lambda x: x.replace(year=year)) + df = df.resample("h").ffill().reindex(timesteps).ffill() + df["year"] = year + df = df.set_index(["year", df.index]) # df.index is timestep + dfs.append(df) + df = pd.concat(dfs) + return df.reindex(n.snapshots) + + +# --------------------------------------------------------------------------- +# Shared construction helpers +# --------------------------------------------------------------------------- + + +def _get_regions_2_add(n: pypsa.Network, flowgates: pd.DataFrame, zone_col: str) -> list[str]: + """Gets regions to add import and export buses to.""" + unique_regions = set(flowgates.r.unique()) | set(flowgates.rr.unique()) + return [x for x in unique_regions if x not in n.buses[zone_col].unique()] + + +def _add_import_export_carriers(n: pypsa.Network, direction: str, co2_emissions: float | None = None) -> None: + """Adds import and export carriers to the network.""" + if direction == "imports": + co2_emissions = 0 if not co2_emissions else co2_emissions + n.add("Carrier", "imports", co2_emissions=co2_emissions, nice_name="Imports") + elif direction == "exports": + n.add("Carrier", "exports", co2_emissions=0, nice_name="Exports") + else: + raise ValueError(f"direction must be either imports or exports; received: {direction}") + + +def _add_import_export_buses(n: pypsa.Network, regions_2_add: list[str], direction: str) -> None: + """Adds import and export buses to the network.""" + if direction == "imports": + suffix = "_imports" + carrier = "imports" + elif direction == "exports": + suffix = "_exports" + carrier = "exports" + else: + raise ValueError(f"direction must be either imports or exports; received: {direction}") + + # cant add in the reeds_state, reeds_zone, reeds_ba, interconnect, trans_reg, trans_grp + # because this information has already been filtered out of the network + + n.add( + "Bus", + regions_2_add, + suffix=suffix, + carrier=carrier, + country=regions_2_add, + ) + + +def _add_import_export_stores(n: pypsa.Network, regions_2_add: list[str], direction: str) -> None: + """Adds import and export stores to the network.""" + if direction == "imports": + n.add( + "Store", + regions_2_add, + bus=[f"{x}_imports" for x in regions_2_add], + suffix="_imports", + carrier="imports", + e_nom=0, + e_nom_extendable=True, + capital_cost=0, + e_nom_min=0, + e_nom_max=1e9, + e_min_pu=-1, + e_max_pu=0, + e_cyclic_per_period=False, + marginal_cost=0, + ) + elif direction == "exports": + n.add( + "Store", + regions_2_add, + bus=[f"{x}_exports" for x in regions_2_add], + suffix="_exports", + carrier="exports", + e_nom_extendable=True, + marginal_cost=0, + e_nom=0, + e_nom_max=1e9, + e_min=0, + e_min_pu=0, + e_max_pu=1, + ) + else: + raise ValueError(f"direction must be either imports or exports; received: {direction}") + + +def _add_import_export_links( + n: pypsa.Network, + flowgates: pd.DataFrame, + fuel_costs: pd.DataFrame | float | str, + direction: str, + zone_col: str = "reeds_zone", + priced: bool = True, +) -> None: + """Adds import and export links to the network. + + ``priced`` is what distinguishes the two representations on the import side: + in ``store`` mode the link carries the import price, in ``generator`` mode + the price sits on the external generator instead and the link is a pure + transfer capacity. + """ + costs = {} + zones_in_model = n.buses[zone_col].dropna().unique() + + for _, row in flowgates.iterrows(): + zone_inside = row.r if row.r in zones_in_model else row.rr + zone_outside = row.r if row.r not in zones_in_model else row.rr + + # extremely crude caching for generating cost timeseries :| + # keyed by the INSIDE zone — checking the outside zone here skipped + # the write whenever an earlier row's inside zone happened to match, + # leaving costs[zone_inside] unset (KeyError on county networks). + if zone_inside not in costs: + if isinstance(fuel_costs, float | int): + costs[zone_inside] = fuel_costs + elif isinstance(fuel_costs, pd.DataFrame): + costs[zone_inside] = _build_cost_timeseries(n, fuel_costs, zone_inside) + else: + costs[zone_inside] = 0 + + marginal_cost = costs[zone_inside] + + capacity = row.value + + """Structre of flowgates is given by: + + r rr value + 0 p6 p8 488.117 + 1 p8 p6 378.458 + 2 p6 p9 4800.000 + ... + """ + + if direction == "imports": + if row.r == zone_inside: # originating at r is exports (ie r -> rr) + continue + name = f"{zone_inside}_{zone_outside}_imports" + bus0 = f"{zone_outside}_imports" + bus1 = zone_inside + carrier = "imports" + if not priced: + marginal_cost = 0 + else: + if row.r == zone_outside: # originating at rr is exports (ie rr -> r) + continue + name = f"{zone_inside}_{zone_outside}_exports" + bus0 = zone_inside + bus1 = f"{zone_outside}_exports" + carrier = "exports" + if isinstance(marginal_cost, pd.Series): + marginal_cost = marginal_cost.mul(-1) # constraint will limit exports + + mc = marginal_cost.value if isinstance(marginal_cost, pd.DataFrame) else marginal_cost + + n.add( + "Link", + name, + bus0=bus0, + bus1=bus1, + carrier=carrier, + p_nom_extendable=False, + p_min_pu=0, + p_max_pu=1, + marginal_cost=mc, + p_nom=capacity, + ) + + +# --------------------------------------------------------------------------- +# `store` representation +# --------------------------------------------------------------------------- + + +def add_elec_imports_exports( + n: pypsa.Network, + direction: str, + flowgates: pd.DataFrame, + fuel_costs: pd.DataFrame | float, + co2_emissions: float = 0, + zone_col: str = "reeds_zone", +): + """Add electricity imports and exports to the network. + + These are capacity constrianed links to/from states outside the model spatial scope. + """ + assert direction in ["imports", "exports"], f"direction must be either imports or exports; received: {direction}" + + regions_2_add = _get_regions_2_add(n, flowgates, zone_col) + _add_import_export_carriers(n, direction, co2_emissions) + _add_import_export_buses(n, regions_2_add, direction) + _add_import_export_stores(n, regions_2_add, direction) + _add_import_export_links(n, flowgates, fuel_costs, direction, zone_col) + + +# --------------------------------------------------------------------------- +# `generator` representation +# --------------------------------------------------------------------------- + + +def inbound_capacity_by_zone( + n: pypsa.Network, + flowgates: pd.DataFrame, + zone_col: str = "reeds_zone", +) -> pd.Series: + """Total interface capacity flowing INTO the footprint, per external zone.""" + zones_in_model = n.buses[zone_col].dropna().unique() + inbound = flowgates[~flowgates.r.isin(zones_in_model) & flowgates.rr.isin(zones_in_model)] + return inbound.groupby("r")["value"].sum() + + +def _internal_zones_served( + n: pypsa.Network, + flowgates: pd.DataFrame, + zone_col: str, +) -> dict[str, list[str]]: + """Map each external zone onto the internal zones it can deliver into.""" + zones_in_model = n.buses[zone_col].dropna().unique() + inbound = flowgates[~flowgates.r.isin(zones_in_model) & flowgates.rr.isin(zones_in_model)] + return {zone: sorted(set(rows.rr)) for zone, rows in inbound.groupby("r")} + + +def _external_generator_cost( + n: pypsa.Network, + fuel_costs: pd.DataFrame | float, + internal_zones: list[str], +) -> pd.Series | float: + """Price for the generic import generator of one external zone. + + The wholesale price table is keyed by the zones INSIDE the model (it is built + from bus attributes), so an external zone has no price of its own. Its + generator is priced at the mean of the internal zones it can reach, which is + the same set of prices the ``store``-mode links would have carried. + """ + if isinstance(fuel_costs, float | int): + return fuel_costs + if not isinstance(fuel_costs, pd.DataFrame): + return 0 + series = [_build_cost_timeseries(n, fuel_costs, zone)["value"] for zone in internal_zones] + series = [s for s in series if not s.isna().all()] + if not series: + return 0 + return pd.concat(series, axis=1).mean(axis=1) + + +def _add_generic_import_generators( + n: pypsa.Network, + regions_2_add: list[str], + flowgates: pd.DataFrame, + fuel_costs: pd.DataFrame | float, + co2_emissions: float, + zone_col: str, +) -> None: + """One priced, non-extendable import generator per external bus.""" + n.add("Carrier", GENERIC_IMPORT_CARRIER, co2_emissions=co2_emissions, nice_name="Unspecified Imports") + + capacity = inbound_capacity_by_zone(n, flowgates, zone_col) + served = _internal_zones_served(n, flowgates, zone_col) + + for zone in regions_2_add: + p_nom = float(capacity.get(zone, 0.0)) + if p_nom <= 0: + logger.info(f"External zone '{zone}' has no inbound interface capacity; no import generator added.") + continue + marginal_cost = _external_generator_cost(n, fuel_costs, served.get(zone, [])) + n.add( + "Generator", + f"{zone}_imports {GENERIC_IMPORT_CARRIER}", + bus=f"{zone}_imports", + carrier=GENERIC_IMPORT_CARRIER, + p_nom=p_nom, + p_nom_extendable=False, + efficiency=1, + marginal_cost=marginal_cost, + ) + + +def map_remote_units_to_zones( + unit_states: pd.Series, + external_zones: list[str], + inbound_capacity: pd.Series, + zone_col: str, + membership: pd.DataFrame | None = None, +) -> pd.Series: + """Assign each contracted remote unit to the external zone it sits behind. + + Candidate zones are the boundary zones whose state matches the unit's + physical state (from ``powerplants.csv``); the candidate with the largest + inbound interface capacity wins. A unit whose state has no direct interface + with the footprint (Utah's Intermountain, say) falls back to the boundary + zone with the largest inbound capacity overall, with a warning. + """ + ranked = [z for z in external_zones if inbound_capacity.get(z, 0.0) > 0] + ranked = sorted(ranked, key=lambda z: (-float(inbound_capacity.get(z, 0.0)), z)) + if not ranked: + raise ValueError("No external zone with inbound interface capacity; cannot place remote contracted units.") + + zone_state = _external_zone_states(ranked, zone_col, membership) + fallback = ranked[0] + + assignment = {} + for unit, state in unit_states.items(): + candidates = [z for z in ranked if zone_state.get(z) == state] + if candidates: + assignment[unit] = candidates[0] + else: + assignment[unit] = fallback + logger.warning( + f"Remote contracted unit '{unit}' sits in state '{state}', which has no direct interface with the " + f"model footprint; placing it behind external zone '{fallback}' (largest inbound capacity).", + ) + return pd.Series(assignment, dtype=object) + + +# County zones are "p" + 5-digit county FIPS; the first two digits are the state. +_STATE_BY_FIPS = { + "01": "AL", + "02": "AK", + "04": "AZ", + "05": "AR", + "06": "CA", + "08": "CO", + "09": "CT", + "10": "DE", + "11": "DC", + "12": "FL", + "13": "GA", + "15": "HI", + "16": "ID", + "17": "IL", + "18": "IN", + "19": "IA", + "20": "KS", + "21": "KY", + "22": "LA", + "23": "ME", + "24": "MD", + "25": "MA", + "26": "MI", + "27": "MN", + "28": "MS", + "29": "MO", + "30": "MT", + "31": "NE", + "32": "NV", + "33": "NH", + "34": "NJ", + "35": "NM", + "36": "NY", + "37": "NC", + "38": "ND", + "39": "OH", + "40": "OK", + "41": "OR", + "42": "PA", + "44": "RI", + "45": "SC", + "46": "SD", + "47": "TN", + "48": "TX", + "49": "UT", + "50": "VT", + "51": "VA", + "53": "WA", + "54": "WV", + "55": "WI", + "56": "WY", +} + +_COUNTY_ZONE_RE = re.compile(r"^p?(\d{2})\d{3}$") + + +def _external_zone_states( + external_zones: list[str], + zone_col: str, + membership: pd.DataFrame | None, +) -> dict[str, str]: + """State code of each external zone. + + With ``topological_boundaries: state`` the zone IS the state code; county + zones ("p" + county FIPS) resolve through the state FIPS prefix; ReEDS + balancing areas resolve through the membership table. + """ + if zone_col == "reeds_state": + return {zone: zone for zone in external_zones} + mbshp = membership.set_index("ba")["st"] if membership is not None else pd.Series(dtype=object) + + def state_of(zone: str) -> str | None: + m = _COUNTY_ZONE_RE.match(zone) + if m and zone not in mbshp.index: + return _STATE_BY_FIPS.get(m.group(1)) + return mbshp.get(zone) + + states = {zone: state_of(zone) for zone in external_zones} + if membership is None and not all(states.values()): + logger.warning("No ReEDS membership table supplied; cannot map external zones onto states.") + return states + + +def attach_remote_contracted_units_externally( + n: pypsa.Network, + bundle: dict, + flowgates: pd.DataFrame, + zone_col: str, + membership: pd.DataFrame | None, +) -> pd.Series: + """Attach the serialized CPUC contracted units at their external buses. + + ``bundle`` is the ``add_electricity`` output described in + :func:`add_electricity.build_remote_unit_bundle`: the fully-derived unit + table (still keyed to the California bus whose VRE profile it borrowed) plus + the borrowed profiles and the cost table needed to resolve capital costs. + Only the ``bus`` column is rewritten here. + """ + units = bundle["units"].copy() + if units.empty: + return pd.Series(dtype=object) + + external_zones = [b[: -len("_imports")] for b in n.buses.index[n.buses.index.str.endswith("_imports")]] + inbound = inbound_capacity_by_zone(n, flowgates, zone_col) + zones = map_remote_units_to_zones(units["state"], external_zones, inbound, zone_col, membership) + + units["zone"] = zones + units["bus"] = zones.map(lambda z: f"{z}_imports") + + for name, unit in units.iterrows(): + logger.info( + f"Remote contracted unit '{name}' ({unit['carrier']}, {unit['p_nom']:.1f} MW, state {unit['state']}) " + f"attached behind external zone '{unit['zone']}' at bus '{unit['bus']}'.", + ) + + add_missing_carriers(n, sorted(set(units["carrier"]))) + attach_remote_units( + n, + units.drop(columns=["zone"]), + bundle["costs"], + bundle["conventional_carriers"], + bundle["unit_commitment"], + profiles=bundle.get("vre_profiles"), + ) + return units["zone"] + + +def _add_generator_representation( + n: pypsa.Network, + flowgates: pd.DataFrame, + fuel_costs: pd.DataFrame | float, + co2_emissions: float, + zone_col: str, + remote_bundle: dict | None, + membership: pd.DataFrame | None, +) -> None: + """Import side of the ``generator`` representation (exports are unchanged).""" + regions_2_add = _get_regions_2_add(n, flowgates, zone_col) + + # emissions move onto the generic import generator's carrier; the `imports` + # carrier now only labels links, which PyPSA never charges emissions to. + _add_import_export_carriers(n, "imports", 0) + _add_import_export_buses(n, regions_2_add, "imports") + _add_generic_import_generators(n, regions_2_add, flowgates, fuel_costs, co2_emissions, zone_col) + _add_import_export_links(n, flowgates, fuel_costs, "imports", zone_col, priced=False) + + if remote_bundle: + attach_remote_contracted_units_externally(n, remote_bundle, flowgates, zone_col, membership) + + +# --------------------------------------------------------------------------- +# Public entry point +# --------------------------------------------------------------------------- + + +def add_external_regions( + n: pypsa.Network, + direction: str, + representation: str, + flowgates: pd.DataFrame, + fuel_costs: pd.DataFrame | float, + co2_emissions: float = 0, + zone_col: str = "reeds_zone", + remote_bundle: dict | None = None, + membership: pd.DataFrame | None = None, +) -> None: + """Add the external-region representation for one direction of trade. + + Parameters + ---------- + direction + ``imports`` or ``exports``. + representation + ``store`` or ``generator`` — see the module docstring. + flowgates + Formatted NARIS flowgates (columns ``r``, ``rr``, ``value``), already + restricted to interfaces that cross the model boundary. + fuel_costs + Output of :func:`resolve_trade_costs`. + remote_bundle, membership + Only used for ``direction='imports'`` in ``generator`` mode: the CPUC + contracted-unit bundle written by ``add_electricity`` and the ReEDS + ``membership.csv`` used to map external zones onto states. + """ + if direction not in ("imports", "exports"): + raise ValueError(f"direction must be either imports or exports; received: {direction}") + if representation not in REPRESENTATIONS: + raise ValueError(f"representation must be one of {REPRESENTATIONS}; received: {representation}") + + if representation == "store" or direction == "exports": + add_elec_imports_exports(n, direction, flowgates, fuel_costs, co2_emissions, zone_col) + return + + _add_generator_representation(n, flowgates, fuel_costs, co2_emissions, zone_col, remote_bundle, membership) + + +def load_remote_unit_bundle(path: str | None) -> dict | None: + """Read the contracted-unit bundle written by ``add_electricity``. + + The rule always declares the file, so a run with contracted resources + disabled (or in ``store`` mode) writes a sentinel ``None``; that is not an + error, it just means there is nothing to place behind the boundary. + """ + if not path: + return None + with open(path, "rb") as f: + bundle = dill.load(f) + if not bundle or bundle.get("units") is None or bundle["units"].empty: + return None + return bundle + + +__all__ = [ + "add_elec_imports_exports", + "add_external_regions", + "calc_import_export_costs", + "convert_flowgates_to_state", + "format_flowgates_for_imports_exports", + "format_import_export_costs", + "inbound_capacity_by_zone", + "load_import_export_costs", + "load_remote_unit_bundle", + "map_remote_units_to_zones", + "resolve_trade_costs", +] diff --git a/workflow/scripts/opts/interfaces.py b/workflow/scripts/opts/interfaces.py index 4dd59c2f..15fcacc1 100644 --- a/workflow/scripts/opts/interfaces.py +++ b/workflow/scripts/opts/interfaces.py @@ -10,7 +10,7 @@ CAISO_Imports,"p9, p10, p11","p2, p5, p6, ...",9728,10208,RESOLVE The caps are applied to the import/export ``Link`` components created by -``add_extra_components.add_elec_imports_exports`` and are therefore a no-op +``external_regions.add_external_regions`` and are therefore a no-op when ``electricity.imports``/``electricity.exports`` are disabled. """ diff --git a/workflow/scripts/plot_statistics.py b/workflow/scripts/plot_statistics.py index 8041386f..54314a24 100644 --- a/workflow/scripts/plot_statistics.py +++ b/workflow/scripts/plot_statistics.py @@ -229,8 +229,8 @@ def plot_capacity_additions_bar( optimal_capacity = optimal_capacity.fillna(0) # Drop the synthetic "imports" carrier — it represents external power - # injection (from trim_network or the imports/exports config), not built - # capacity, and would otherwise dominate the bar by orders of magnitude. + # injection (from the imports/exports config), not built capacity, and would + # otherwise dominate the bar by orders of magnitude. hidden_carriers = {"imports", "Imports"} optimal_capacity = optimal_capacity.drop( index=[c for c in hidden_carriers if c in optimal_capacity.index], diff --git a/workflow/scripts/test/test_interfaces.py b/workflow/scripts/test/test_interfaces.py index 91b0f9c5..6035afbd 100644 --- a/workflow/scripts/test/test_interfaces.py +++ b/workflow/scripts/test/test_interfaces.py @@ -1,13 +1,16 @@ """ -Test the aggregate transmission interface limits. +Test the aggregate transmission interface limits and the external-region build. This module contains tests for the RESOLVE/NARIS style interface constraints -applied to the electricity import/export links in PyPSA-USA. +applied to the electricity import/export links in PyPSA-USA, and for the two +representations of out-of-footprint supply built by ``external_regions``. """ +import logging import os import sys +import numpy as np import pandas as pd import pypsa import pytest @@ -15,6 +18,12 @@ sys.path.append(os.path.join(os.path.dirname(__file__), "..")) from _helpers import get_multiindex_snapshots +from external_regions import ( + GENERIC_IMPORT_CARRIER, + add_external_regions, + inbound_capacity_by_zone, + map_remote_units_to_zones, +) from opts.interfaces import ( _boundary_links, _parse_regions, @@ -36,7 +45,7 @@ def interface_network(): """ Build a small network with electricity import and export links. - Mirrors the conventions of ``add_extra_components.add_elec_imports_exports``: + Mirrors the conventions of ``external_regions.add_elec_imports_exports``: external buses are named ``{zone}_imports`` / ``{zone}_exports`` with the matching carrier, import links run from the external bus into the model and export links run the other way. @@ -285,3 +294,349 @@ def extra_functionality(n, sns): n.optimize(solver_name="glpk", multi_investment_periods=True, extra_functionality=extra_functionality) assert not [c for c in n.model.constraints if c.startswith("interface_limit-")] + + +# --------------------------------------------------------------------------- +# external_regions: the two import representations +# --------------------------------------------------------------------------- + +MEMBERSHIP = pd.DataFrame( + [ + {"ba": "p9", "st": "CA"}, + {"ba": "p12", "st": "NV"}, + {"ba": "p13", "st": "NV"}, + {"ba": "p28", "st": "AZ"}, + ], +) + +# Mirrors the CA boundary: two NV zones (p13 much larger than p12), one AZ zone, +# and one outbound-only interface so the inbound/outbound split is exercised. +FLOWGATES = pd.DataFrame( + [ + {"r": "p13", "rr": "p9", "value": 2000.0}, + {"r": "p12", "rr": "p9", "value": 500.0}, + {"r": "p28", "rr": "p9", "value": 1000.0}, + {"r": "p9", "rr": "p13", "value": 1500.0}, + ], +) + + +@pytest.fixture +def footprint_network(): + """A one-zone footprint (California's p9) with nothing external attached yet.""" + n = pypsa.Network() + n.snapshots = get_multiindex_snapshots( + sns_config={"start": "2030-01-01 00:00", "end": "2030-01-01 03:00", "inclusive": "both"}, + invest_periods=[2030], + ) + n.set_investment_periods(periods=[2030]) + + n.add("Carrier", "AC", co2_emissions=0) + n.add("Carrier", "solar", co2_emissions=0) + n.add( + "Bus", + "p9", + carrier="AC", + country="p9", + interconnect="western", + reeds_state="CA", + reeds_zone="p9", + ) + n.add( + "Generator", + "p9 solar", + bus="p9", + carrier="solar", + p_nom=100, + p_max_pu=pd.Series(0.5, index=n.snapshots), + ) + return n + + +def wholesale_costs(): + """A `format_import_export_costs`-shaped price table for the internal zone.""" + return pd.DataFrame( + {"zone": ["p9"], "value": [42.0], "units": ["usd/mwh"]}, + index=pd.to_datetime(["2030-01-01"]), + ) + + +def import_links(n): + return n.links[n.links.carrier == "imports"] + + +def test_inbound_capacity_by_zone_ignores_outbound_rows(footprint_network): + inbound = inbound_capacity_by_zone(footprint_network, FLOWGATES, "reeds_zone") + assert inbound.to_dict() == {"p12": 500.0, "p13": 2000.0, "p28": 1000.0} + + +def test_generator_mode_builds_external_buses_generators_and_links(footprint_network): + n = footprint_network + add_external_regions(n, "imports", "generator", FLOWGATES, 30.0, co2_emissions=0.428, zone_col="reeds_zone") + + for zone in ("p12", "p13", "p28"): + assert f"{zone}_imports" in n.buses.index + assert n.buses.at[f"{zone}_imports", "carrier"] == "imports" + + gens = n.generators[n.generators.carrier == GENERIC_IMPORT_CARRIER] + assert sorted(gens.bus) == ["p12_imports", "p13_imports", "p28_imports"] + assert not gens.p_nom_extendable.any() + # p_nom is the zone's total INBOUND interface capacity + assert gens.set_index("bus")["p_nom"].to_dict() == { + "p12_imports": 500.0, + "p13_imports": 2000.0, + "p28_imports": 1000.0, + } + assert (gens.marginal_cost == 30.0).all() + + links = import_links(n) + assert sorted(links.index) == ["p9_p12_imports", "p9_p13_imports", "p9_p28_imports"] + assert links.set_index("bus0")["p_nom"].to_dict() == { + "p12_imports": 500.0, + "p13_imports": 2000.0, + "p28_imports": 1000.0, + } + assert (links.bus1 == "p9").all() + assert not links.p_nom_extendable.any() + + # generator mode prices the GENERATOR, so the interface link is free + assert (links.marginal_cost == 0).all() + # and there is no bottomless import Store + assert n.stores[n.stores.carrier == "imports"].empty + + +def test_store_mode_prices_the_links_and_keeps_the_store(footprint_network): + n = footprint_network + add_external_regions(n, "imports", "store", FLOWGATES, 30.0, co2_emissions=0.428, zone_col="reeds_zone") + + links = import_links(n) + assert (links.marginal_cost == 30.0).all() + assert sorted(n.stores[n.stores.carrier == "imports"].index) == [ + "p12_imports", + "p13_imports", + "p28_imports", + ] + assert n.generators[n.generators.carrier == GENERIC_IMPORT_CARRIER].empty + + +def test_representations_agree_on_interface_capacity(footprint_network): + """The transfer capacity a zone can deliver is the same in both modes.""" + store = footprint_network.copy() + add_external_regions(store, "imports", "store", FLOWGATES, 30.0, zone_col="reeds_zone") + generator = footprint_network.copy() + add_external_regions(generator, "imports", "generator", FLOWGATES, 30.0, zone_col="reeds_zone") + + pd.testing.assert_series_equal( + import_links(store).set_index("bus0")["p_nom"].sort_index(), + import_links(generator).set_index("bus0")["p_nom"].sort_index(), + ) + + +def test_generator_mode_moves_emissions_onto_the_generator_carrier(footprint_network): + n = footprint_network + add_external_regions(n, "imports", "generator", FLOWGATES, 30.0, co2_emissions=0.428, zone_col="reeds_zone") + + # links never carry emissions in PyPSA, so the factor must sit on the generator + assert n.carriers.at[GENERIC_IMPORT_CARRIER, "co2_emissions"] == pytest.approx(0.428) + assert n.carriers.at["imports", "co2_emissions"] == 0 + + +def test_store_mode_keeps_emissions_on_the_imports_carrier(footprint_network): + n = footprint_network + add_external_regions(n, "imports", "store", FLOWGATES, 30.0, co2_emissions=0.428, zone_col="reeds_zone") + assert n.carriers.at["imports", "co2_emissions"] == pytest.approx(0.428) + + +def test_generator_price_comes_from_the_wholesale_table(footprint_network): + n = footprint_network + add_external_regions(n, "imports", "generator", FLOWGATES, wholesale_costs(), zone_col="reeds_zone") + + name = f"p13_imports {GENERIC_IMPORT_CARRIER}" + assert name in n.generators_t.marginal_cost.columns + assert n.generators_t.marginal_cost[name].to_numpy() == pytest.approx(42.0) + + +@pytest.mark.parametrize("representation", ["store", "generator"]) +def test_exports_are_identical_across_representations(footprint_network, representation): + """Export construction is deliberately untouched by the representation switch.""" + n = footprint_network + add_external_regions(n, "exports", representation, FLOWGATES, -30.0, zone_col="reeds_zone") + + links = n.links[n.links.carrier == "exports"] + assert sorted(links.index) == ["p9_p13_exports"] + assert links.at["p9_p13_exports", "bus0"] == "p9" + assert links.at["p9_p13_exports", "bus1"] == "p13_exports" + assert links.at["p9_p13_exports", "p_nom"] == 1500.0 + # the export revenue stays on the link in both modes + assert links.at["p9_p13_exports", "marginal_cost"] == -30.0 + assert not n.stores[n.stores.carrier == "exports"].empty + # nothing but export links may inject into the export bus + assert n.generators[n.generators.bus == "p13_exports"].empty + + +# --------------------------------------------------------------------------- +# external_regions: contracted units behind the boundary +# --------------------------------------------------------------------------- + + +def test_map_remote_units_to_zones_prefers_the_largest_matching_zone(): + inbound = pd.Series({"p12": 500.0, "p13": 2000.0, "p28": 1000.0}) + zones = map_remote_units_to_zones( + pd.Series({"R Hoover": "NV", "R Palo Verde": "AZ"}), + ["p12", "p13", "p28"], + inbound, + "reeds_zone", + MEMBERSHIP, + ) + assert zones["R Hoover"] == "p13" # NV, and p13 > p12 + assert zones["R Palo Verde"] == "p28" + + +def test_map_remote_units_to_zones_falls_back_and_warns(caplog): + inbound = pd.Series({"p12": 500.0, "p13": 2000.0, "p28": 1000.0}) + with caplog.at_level(logging.WARNING): + zones = map_remote_units_to_zones( + pd.Series({"R Intermountain": "UT"}), + ["p12", "p13", "p28"], + inbound, + "reeds_zone", + MEMBERSHIP, + ) + assert zones["R Intermountain"] == "p13" # largest inbound capacity overall + assert "R Intermountain" in caplog.text + assert "UT" in caplog.text + + +def test_map_remote_units_to_zones_resolves_county_fips_states(): + # County zones are "p" + county FIPS and are absent from the BA membership + # table; the state must come from the FIPS prefix (32=NV, 04=AZ). + inbound = pd.Series({"p32003": 3000.0, "p04012": 1000.0, "p04027": 2000.0}) + zones = map_remote_units_to_zones( + pd.Series({"R Hoover": "NV", "R Palo Verde": "AZ"}), + ["p32003", "p04012", "p04027"], + inbound, + "county", + MEMBERSHIP, + ) + assert zones["R Hoover"] == "p32003" + assert zones["R Palo Verde"] == "p04027" # AZ, and p04027 > p04012 + + +def test_map_remote_units_to_zones_state_boundaries_need_no_membership(): + inbound = pd.Series({"NV": 2000.0, "AZ": 1000.0}) + zones = map_remote_units_to_zones( + pd.Series({"R Apex": "NV"}), + ["NV", "AZ"], + inbound, + "reeds_state", + None, + ) + assert zones["R Apex"] == "NV" + + +def remote_bundle(n): + """A toy `build_remote_unit_bundle` output: one AZ firm unit, one UT VRE unit.""" + units = pd.DataFrame( + [ + { + "name": "R Palo Verde", + "carrier": "nuclear", + "bus": "p9", + "p_nom": 600.0, + "state": "AZ", + "efficiency": 0.33, + "marginal_cost": 8.0, + "heat_rate": 10.0, + "summer_derate": 1.0, + "winter_derate": 1.0, + "ramp_limit_up": 1.0, + "ramp_limit_down": 1.0, + "min_up_time": 0, + "min_down_time": 0, + "start_up_cost": 0.0, + "fuel_cost": 1.0, + "min_load_pu": 0.0, + "build_year": 1988, + "duration": 4.0, + }, + { + "name": "R Cape Solar", + "carrier": "solar", + "bus": "p9", + "p_nom": 50.0, + "state": "UT", + "efficiency": 1.0, + "marginal_cost": 0.0, + "heat_rate": 0.0, + "summer_derate": 1.0, + "winter_derate": 1.0, + "ramp_limit_up": 1.0, + "ramp_limit_down": 1.0, + "min_up_time": 0, + "min_down_time": 0, + "start_up_cost": 0.0, + "fuel_cost": 0.0, + "min_load_pu": 0.0, + "build_year": 2024, + "duration": 4.0, + }, + ], + ).set_index("name") + return { + "units": units, + "vre_profiles": pd.DataFrame({"R Cape Solar": np.full(len(n.snapshots), 0.4)}, index=n.snapshots), + "costs": pd.DataFrame(), + "conventional_carriers": ["nuclear"], + "unit_commitment": False, + } + + +def test_remote_bundle_attaches_behind_the_matching_external_bus(footprint_network, caplog): + n = footprint_network + with caplog.at_level(logging.WARNING): + add_external_regions( + n, + "imports", + "generator", + FLOWGATES, + 30.0, + co2_emissions=0.428, + zone_col="reeds_zone", + remote_bundle=remote_bundle(n), + membership=MEMBERSHIP, + ) + + # AZ unit lands behind the AZ boundary zone + assert n.generators.at["R Palo Verde", "bus"] == "p28_imports" + assert n.generators.at["R Palo Verde", "p_nom"] == 600.0 + assert not n.generators.at["R Palo Verde", "p_nom_extendable"] + + # UT has no direct CA interface -> fallback to the largest inbound zone + assert n.generators.at["R Cape Solar", "bus"] == "p13_imports" + assert "R Cape Solar" in caplog.text + + # the borrowed profile travels with the bundle + assert n.generators_t.p_max_pu["R Cape Solar"].to_numpy() == pytest.approx(0.4) + + # and none of them sit inside the footprint any more + assert "p9" not in set(n.generators.loc[["R Palo Verde", "R Cape Solar"], "bus"]) + + +def test_remote_bundle_is_ignored_in_store_mode(footprint_network): + n = footprint_network + add_external_regions( + n, + "imports", + "store", + FLOWGATES, + 30.0, + zone_col="reeds_zone", + remote_bundle=remote_bundle(n), + membership=MEMBERSHIP, + ) + assert "R Palo Verde" not in n.generators.index + + +def test_unknown_representation_raises(footprint_network): + with pytest.raises(ValueError, match="representation"): + add_external_regions(footprint_network, "imports", "banana", FLOWGATES, 30.0, zone_col="reeds_zone")