fix(rps): credit behind-the-meter rooftop solar against RPS targets

docs/source/data-sectors.md

@@ -89,7 +89,7 @@ Within the transport sector, demand response is applied to the aggregation bus d
 
 The diagram below illustrates the workflow of PyPSA-USA Sector. Many rules overlap with the electricity sector workflow; however, several additional rules are also present.
 
-:::{figure-md} workflow
+:::{figure-md} sector-workflow
 <img src="./_static/dag_sector.jpg" width="700px">
 
 Snakemake DAG for sector coupled studies

docs/source/data-transmission.md

@@ -34,7 +34,7 @@ While representative of the US electricity system, the TAMU network is synthetic
 
 
 
-```{info}
+```{note}
 See the [Spatial Configuration](./config-spatial.md) page for information on how to choose between networks.
 ```
 

docs/source/publications.bib

@@ -1,3 +1,33 @@
+@misc{barnes2026b,
+  title = {Near-Term Emission Targets Need Immediate Attention in the USA},
+  author = {Trevor Barnes and Kamran Tehranchi and Brad Reinholz and Malcolm Metcalfe and Taco Niet},
+  year = {2026},
+  archivePrefix = {arXiv},
+  url = {https://arxiv.org/abs/2607.01471},
+  doi = {10.48550/arXiv.2607.01471},
+}
+
+@misc{ai2026,
+  title = {Economic Valuation and Optimal Deployment of Static Synchronous Series Compensators for U.S. Power System Expansion},
+  author = {Wei Ai and Vladimir Dvorkin and Michael T. Craig},
+  year = {2026},
+  archivePrefix = {arXiv},
+  url = {https://arxiv.org/abs/2605.00734},
+  doi = {10.48550/arXiv.2605.00734},
+}
+
+@article{barnes_2026a,
+    doi = {10.1371/journal.pclm.0000918},
+    url = {https://journals.plos.org/climate/article?id=10.1371/journal.pclm.0000918},
+    author = {Barnes, Trevor and Tehranchi, Kamran and Reinholz, Bradley and Metcalfe, Malcolm and Niet, Taco},
+    title = {Multi-sector demand response for cost optimal energy transitions},
+    year = {2026},
+    volume = {5},
+    issue = {5},
+    pages = {e0000918},
+    journal = {PLOS Climate}
+}
+
 @misc{tehranchi_2024,
     doi = {10.2139/ssrn.5029120},
     url = {https://www.ssrn.com/abstract=5029120},

workflow/repo_data/policy_constraints/small_scale_solar.csv

@@ -0,0 +1,49 @@
+state,year,generation_mwh
+AL,2023,460000
+AZ,2023,5890000
+AR,2023,190000
+CA,2023,39800000
+CO,2023,2640000
+CT,2023,1530000
+DE,2023,390000
+FL,2023,9620000
+GA,2023,2180000
+ID,2023,480000
+IL,2023,1130000
+IN,2023,560000
+IA,2023,290000
+KS,2023,270000
+KY,2023,200000
+LA,2023,460000
+ME,2023,820000
+MD,2023,2050000
+MA,2023,4870000
+MI,2023,1800000
+MN,2023,1010000
+MS,2023,150000
+MO,2023,470000
+MT,2023,170000
+NE,2023,200000
+NV,2023,2100000
+NH,2023,680000
+NJ,2023,4280000
+NM,2023,700000
+NY,2023,4550000
+NC,2023,2060000
+ND,2023,50000
+OH,2023,1280000
+OK,2023,340000
+OR,2023,1460000
+PA,2023,1870000
+RI,2023,520000
+SC,2023,680000
+SD,2023,70000
+TN,2023,480000
+TX,2023,8700000
+UT,2023,1090000
+VT,2023,610000
+VA,2023,1740000
+WA,2023,1110000
+WV,2023,80000
+WI,2023,760000
+WY,2023,120000

workflow/rules/retrieve.smk

@@ -232,8 +232,7 @@ rule retrieve_caiso_data:
         fuel_prices=DATA + "costs/caiso_ng_power_prices.csv",
     log:
         LOGS + "retrieve_caiso_data.log",
-    shadow:
-        "minimal"
+    # shadow: "minimal"  # disabled on Windows (symlink creation requires Developer Mode)
     resources:
         walltime="00:10:00",
         mem_mb=2000,
@@ -254,6 +253,23 @@ rule retrieve_pudl:
         "../scripts/retrieve_pudl.py"
 
 
+rule retrieve_small_scale_solar:
+    params:
+        eia_api=config.get("api", {}).get("eia", None),
+        planning_horizons=config["scenario"]["planning_horizons"],
+    input:
+        fallback="repo_data/policy_constraints/small_scale_solar.csv",
+    output:
+        small_scale_solar=DATA + "eia/small_scale_solar.csv",
+    log:
+        LOGS + "retrieve_small_scale_solar.log",
+    resources:
+        walltime="00:10:00",
+        mem_mb=1000,
+    script:
+        "../scripts/retrieve_small_scale_solar.py"
+
+
 if "EGS" in config["electricity"]["extendable_carriers"]["Generator"]:
 
     rule retrieve_egs:

workflow/rules/solve_electricity.smk

@@ -29,6 +29,7 @@ rule solve_network:
         safer_reeds="config/policy_constraints/reeds/prm_annual.csv",
         rps_reeds="config/policy_constraints/reeds/rps_fraction.csv",
         ces_reeds="config/policy_constraints/reeds/ces_fraction.csv",
+        small_scale_solar=DATA + "eia/small_scale_solar.csv",
         pop_layout=pop_layout_input,
         ev_policy=ev_policy_input,
     output:

workflow/scripts/eia.py

@@ -490,6 +490,43 @@ def data_creator(self):
         return _ElectricPowerOperationalData(self.sector, self.year, self.api_key)
 
 
+class SmallScaleSolar(EiaData):
+    """State-level small-scale (behind-the-meter) solar PV annual generation.
+
+    Fetches sector-98 (below 1 MW nameplate) solar generation from the EIA
+    Electric Power Operational Data API for a range of years.
+
+    Parameters
+    ----------
+    start_year : int
+        First year to retrieve (EIA data available from ~2014).
+    end_year : int
+        Last year to retrieve. Capped at the latest available year.
+    api : str
+        EIA API v2 key.
+
+    Examples
+    --------
+    >>> df = SmallScaleSolar(2014, 2023, "YOUR_API_KEY").get_data()
+    """
+
+    LAST_AVAILABLE_YEAR: ClassVar[int] = 2023
+
+    def __init__(self, start_year: int, end_year: int, api: str) -> None:
+        self.start_year = start_year
+        self.end_year = min(end_year, self.LAST_AVAILABLE_YEAR)
+        self.api = api
+        if end_year > self.LAST_AVAILABLE_YEAR:
+            logger.warning(
+                f"Small-scale solar data only available through {self.LAST_AVAILABLE_YEAR}. "
+                f"Capping end_year at {self.LAST_AVAILABLE_YEAR}.",
+            )
+
+    def data_creator(self):
+        """Initializes data extractor."""
+        return _SmallScaleSolarData(self.start_year, self.end_year, self.api)
+
+
 # product
 class DataExtractor(ABC):
     """Extracts and formats data."""
@@ -1839,6 +1876,52 @@ def format_data(self, df: pd.DataFrame) -> pd.DataFrame:
         return self._assign_dtypes(df)
 
 
+class _SmallScaleSolarData(DataExtractor):
+    """Fetches annual small-scale (behind-the-meter) solar generation by state.
+
+    Uses the EIA Electric Power Operational Data endpoint with sector 98
+    (small-scale photovoltaic, below 1 MW nameplate) and fuel type SUN.
+    Generation values are returned in MWh (converted from the API's thousand MWh).
+    """
+
+    SECTOR_ID: ClassVar[int] = 98
+    FUEL_TYPE: ClassVar[str] = "SUN"
+
+    def __init__(self, start_year: int, end_year: int, api_key: str) -> None:
+        super().__init__(end_year, api_key)  # self.year = end_year
+        self.start_year = start_year
+
+    def build_url(self) -> str:
+        base_url = "electricity/electric-power-operational-data/data/"
+        facets = (
+            f"frequency=annual&data[0]=generation"
+            f"&facets[fueltypeid][]={self.FUEL_TYPE}"
+            f"&facets[sectorid][]={self.SECTOR_ID}"
+            f"&start={self.start_year}&end={self.year}"
+            f"&sort[0][column]=period&sort[0][direction]=asc"
+            f"&offset=0&length=5000"
+        )
+        return f"{API_BASE}{base_url}?api_key={self.api_key}&{facets}"
+
+    def format_data(self, df: pd.DataFrame) -> pd.DataFrame:
+        df = df.rename(
+            columns={
+                "location": "state",
+                "period": "year",
+                "generation": "generation_mwh",
+            },
+        )
+        df = df[["state", "year", "generation_mwh"]].copy()
+        df["year"] = df["year"].astype(int)
+        # EIA reports generation in thousand MWh; convert to MWh
+        df["generation_mwh"] = pd.to_numeric(df["generation_mwh"], errors="coerce") * 1_000
+        df = df.dropna(subset=["generation_mwh"])
+        df = df[df["generation_mwh"] > 0]
+        # Exclude national aggregates (e.g., "US") and sub-state numeric codes
+        df = df[(df["state"].str.len() == 2) & (df["state"] != "US")]
+        return df.sort_values(["state", "year"]).reset_index(drop=True)
+
+
 if __name__ == "__main__":
     with open("./../config/config.api.yaml") as file:
         yaml_data = yaml.safe_load(file)

workflow/scripts/opts/policy.py

@@ -19,6 +19,7 @@
     "offwind",
     "offwind_floating",
     "solar",
+    "solar-rooftop",
     "hydro",
     "geothermal",
     "biomass",
@@ -247,18 +248,54 @@ def _collapse_portfolio_standards(n: pypsa.Network, planning_horizons: list[int]
     return portfolio_standards
 
 
-def add_RPS_constraints(n, config, snakemake=None):
+def _load_small_scale_solar(snakemake) -> pd.Series:
+    """
+    Load state-level small-scale (behind-the-meter) solar generation.
+
+    Returns a Series indexed by (state, year) with generation in MWh, or an
+    empty Series if the input is not available on the snakemake object.
+    """
+    path = getattr(snakemake.input, "small_scale_solar", None)
+    if path is None:
+        return pd.Series(dtype=float)
+
+    df = pd.read_csv(path, dtype={"state": str, "year": int, "generation_mwh": float})
+    return df.set_index(["state", "year"])["generation_mwh"]
+
+
+def add_RPS_constraints(n, config, snakemake=None, sector=False):
     """
     Add Renewable Portfolio Standards (RPS) constraints to the network.
 
     This function enforces constraints on the percentage of electricity generation
     from renewable energy sources for specific regions and planning horizons.
     It reads the necessary data from configuration files and the network.
 
-    The differenct between electrical and sector implementation is:
+    The difference between electrical and sector implementation is:
     - Electrical applies RPS against exogenously defined demand
     - Sector applies RPS against endogenously solved power sector generation
 
+    When ``snakemake.input.small_scale_solar`` is provided, the demand basis for
+    each constraint is adjusted from net load to gross load by adding back the
+    behind-the-meter (rooftop) solar generation that is embedded as a demand
+    reduction in the EIA 930 input data.  This ensures that existing rooftop
+    solar receives credit toward the RPS target even when it is not explicitly
+    modelled as a Generator in the network.
+
+    The adjusted RHS is:
+        rhs = pct * net_load - (1 - pct) * rooftop_gen
+            = pct * gross_load - rooftop_gen
+
+    **Demand source compatibility:** The BTM credit is only appropriate when
+    the network load time-series is derived from EIA 930 *net* generation data
+    (``demand.profile: eia``), in which case behind-the-meter solar is already
+    subtracted from the reported load.  If the demand profile is sourced from
+    EFS or AEO projections — which typically report *gross* electricity sales
+    and do not subtract BTM generation — the BTM credit should not be applied.
+    To disable it, simply omit the ``small_scale_solar`` input from the
+    ``solve_network`` rule (or leave ``api.eia`` unconfigured so the fallback
+    CSV is not forward-filled beyond its data year).
+
     Parameters
     ----------
     n : pypsa.Network
@@ -299,21 +336,58 @@ def add_RPS_constraints(n, config, snakemake=None):
         ces_reeds,
     )
 
+    # Small-scale solar generation by (state, year) in MWh — may be empty if
+    # the input file was not provided.
+    small_scale_solar = _load_small_scale_solar(snakemake)
+    using_btm_credit = not small_scale_solar.empty
+
     for _, constraint_row in portfolio_standards.iterrows():
         region_list = [region.strip() for region in constraint_row.region.split(",")]
         region_buses = get_region_buses(n, region_list)
         if region_buses.empty:
             continue
 
+        # Net load from the grid (EIA 930 data already has BTM solar subtracted)
         region_demand = (
             n.loads_t.p_set.loc[constraint_row.planning_horizon]
             .loc[:, n.loads.bus.isin(region_buses.index)]
             .sum()
             .sum()
         )
-        region_rps_rhs = int(constraint_row.pct * region_demand)
+
+        # Credit existing behind-the-meter rooftop solar toward the RPS target.
+        #
+        # Derivation:
+        #   statutory target: utility_renewables + rooftop >= pct * gross_load
+        #   gross_load = net_load + rooftop
+        #   rearranged:  utility_renewables >= pct * net_load - (1 - pct) * rooftop
+        #
+        # When rooftop solar is explicitly modelled as a Generator (carrier
+        # "solar-rooftop") in the network, it already appears on the LHS and
+        # n.loads_t.p_set is the gross load, so no adjustment is needed.
+        # The BTM credit here only applies to the residual rooftop generation
+        # that is embedded as a demand reduction and has no generator in the
+        # network.
+        rooftop_gen = 0.0
+        if using_btm_credit:
+            for state in region_list:
+                key = (state.strip(), constraint_row.planning_horizon)
+                if key in small_scale_solar.index:
+                    rooftop_gen += small_scale_solar[key]
+
+        pct = constraint_row.pct
+        region_rps_rhs = max(int(pct * region_demand - (1 - pct) * rooftop_gen), 0)
+
         portfolio_standards.loc[constraint_row.name, "rps_rhs"] = region_rps_rhs
 
+        if using_btm_credit and rooftop_gen > 0:
+            logger.info(
+                f"RPS demand basis for {constraint_row.region} ({constraint_row.planning_horizon}): "
+                f"net_load={region_demand / 1e6:.1f} TWh, "
+                f"btm_rooftop={rooftop_gen / 1e6:.2f} TWh, "
+                f"adjusted_rhs={region_rps_rhs / 1e6:.1f} TWh",
+            )
+
     # Iterate through constraints and add RPS constraints to the model
     for (rec_trading_zone, planning_horizon, policy_carriers), zone_constraints in portfolio_standards.groupby(
         ["rec_trading_zone", "planning_horizon", "carrier"],
@@ -328,7 +402,7 @@ def add_RPS_constraints(n, config, snakemake=None):
         region_gens_eligible = region_gens[region_gens.carrier.isin(carriers)]
 
         if region_gens_eligible.empty:
-            return
+            continue  # skip this constraint group; do not exit the whole function
 
         # Eligible generation
         p_eligible = n.model["Generator-p"].sel(