|
| 1 | +import logging |
| 2 | +import os |
| 3 | +import sys |
| 4 | + |
| 5 | +sys.path.append(os.path.abspath(os.path.dirname(__file__))) |
| 6 | +sys.path.append( |
| 7 | + os.path.abspath(os.path.join(os.path.dirname(__file__), "..")) |
| 8 | +) |
| 9 | + |
| 10 | +import matplotlib.pyplot as plt |
| 11 | +import numpy as np |
| 12 | +import pandas as pd |
| 13 | +import pypsa |
| 14 | +import os, re, collections |
| 15 | +import itertools |
| 16 | + |
| 17 | +from _helpers import configure_logging, mock_snakemake |
| 18 | + |
| 19 | +groups = { |
| 20 | + "gas": ["gas CHP", "OCGT", "CCGT", "gas"], |
| 21 | + "heat vent": ["heat vent"], |
| 22 | + "water tanks": ["water tank", "water pit"], |
| 23 | + "heat pump" : ["heat pump"], |
| 24 | + "resistive heater" : ["resistive heater"], |
| 25 | + "biomass": ["biomass"], |
| 26 | + "lignite": ["lignite"], |
| 27 | + "coal": ["coal"], |
| 28 | + "oil": ["oil"], |
| 29 | + "waste": ["waste"], |
| 30 | + "solar": ["solar"], |
| 31 | + "offwind": ["offwind"], |
| 32 | +} |
| 33 | + |
| 34 | +def aggregate_by_keywords(opex_comp_agg, groups): |
| 35 | + """ |
| 36 | + Aggregate rows in opex_comp_agg according to keyword groups. |
| 37 | + |
| 38 | + Parameters |
| 39 | + ---------- |
| 40 | + opex_comp_agg : pd.DataFrame |
| 41 | + DataFrame with row index as technology names. |
| 42 | + groups : dict |
| 43 | + Keys = new aggregated name, |
| 44 | + Values = list of substrings to match in the index. |
| 45 | + |
| 46 | + Returns |
| 47 | + ------- |
| 48 | + pd.DataFrame |
| 49 | + """ |
| 50 | + df_out = opex_comp_agg.copy() |
| 51 | + for new_name, keywords in groups.items(): |
| 52 | + mask = df_out.index.to_series().str.contains("|".join(keywords)) |
| 53 | + if mask.any(): |
| 54 | + summed = df_out.loc[mask].sum() |
| 55 | + df_out = df_out.drop(df_out.index[mask]) |
| 56 | + df_out.loc[new_name] = summed |
| 57 | + return df_out |
| 58 | + |
| 59 | + |
| 60 | + |
| 61 | +if __name__ == "__main__": |
| 62 | + if "snakemake" not in globals(): |
| 63 | + import os |
| 64 | + import sys |
| 65 | + |
| 66 | + from _helpers import mock_snakemake |
| 67 | + |
| 68 | + snakemake = mock_snakemake( |
| 69 | + "regret_plots", |
| 70 | + ) |
| 71 | + |
| 72 | + configure_logging(snakemake) |
| 73 | + config = snakemake.config |
| 74 | + planning_horizons = snakemake.params.planning_horizons |
| 75 | + scenarios = ["AriadneDemand", "LowDemand"] |
| 76 | + tech_colors = snakemake.params.plotting["tech_colors"] |
| 77 | + |
| 78 | + # Nested dict: networks[year][scenario][decision] = Network |
| 79 | + networks = collections.defaultdict( |
| 80 | + lambda: collections.defaultdict(dict) |
| 81 | + ) |
| 82 | + |
| 83 | + for fn in snakemake.input.regret_networks: |
| 84 | + parts = fn.split(os.sep) |
| 85 | + |
| 86 | + # scenario is the folder name 2 levels up |
| 87 | + scenario = parts[-3] |
| 88 | + if scenario not in scenarios: |
| 89 | + raise ValueError(f"Unexpected scenario '{scenario}' in {fn}. Allowed: {scenarios}") |
| 90 | + |
| 91 | + # extract year (4 digits before .nc) |
| 92 | + m = re.search(r"_(\d{4})\.nc$", fn) |
| 93 | + if not m: |
| 94 | + raise ValueError(f"Could not parse year from {fn}") |
| 95 | + year = int(m.group(1)) |
| 96 | + |
| 97 | + # extract decision_* (string until the 2nd underscore in filename) |
| 98 | + filename = parts[-1] |
| 99 | + decision = "_".join(filename.split("_")[:2]) |
| 100 | + if not decision.startswith("decision"): |
| 101 | + raise ValueError(f"Unexpected decision string in {filename}") |
| 102 | + |
| 103 | + # load and store |
| 104 | + # print(f"Loading {fn} ...") |
| 105 | + # print(f" scenario: {scenario}, year: {year}, decision: {decision}") |
| 106 | + networks[year][scenario][decision] = pypsa.Network(fn) |
| 107 | + |
| 108 | + # ensure output directory exist |
| 109 | + for dir in snakemake.output[-1]: |
| 110 | + if not os.path.exists(dir): |
| 111 | + os.makedirs(dir) |
| 112 | + |
| 113 | + # Plot electricity price duration curves |
| 114 | + |
| 115 | + fig, ax = plt.subplots(figsize=(10, 15), nrows=3, ncols=1) |
| 116 | + ax = ax.flatten() |
| 117 | + |
| 118 | + years = [2025, 2030, 2035] |
| 119 | + scenarios = ["AriadneDemand", "LowDemand"] |
| 120 | + decisions = ["decision_AriadneDemand", "decision_LowDemand"] |
| 121 | + |
| 122 | + for i, year in enumerate(years): |
| 123 | + for scenario, decision in itertools.product(scenarios, decisions): |
| 124 | + |
| 125 | + n = networks[year][scenario][decision] |
| 126 | + lmps = n.buses_t.marginal_price.loc[:, |
| 127 | + (n.buses.carrier == "AC") & (n.buses.index.str.startswith("DE"))] |
| 128 | + lmps_sorted = pd.DataFrame(lmps.values.flatten(), columns=["lmp"]).sort_values(by="lmp", ascending=False) |
| 129 | + lmps_sorted["percentage"] = np.arange(len(lmps_sorted)) / len(lmps_sorted) * 100 |
| 130 | + |
| 131 | + ax[i].plot( |
| 132 | + lmps_sorted["percentage"], |
| 133 | + lmps_sorted["lmp"], |
| 134 | + label=f"{scenario}_{decision} (avg: {lmps_sorted['lmp'].mean():.2f})" |
| 135 | + ) |
| 136 | + |
| 137 | + ax[i].legend() |
| 138 | + ax[i].set_xlabel("Percentage of time") |
| 139 | + ax[i].set_ylabel("€/MWh") |
| 140 | + ax[i].set_title(f"Price duration curves {year}") |
| 141 | + |
| 142 | + plt.tight_layout() |
| 143 | + plt.savefig(snakemake.output.elec_price_comp_de, bbox_inches="tight") |
| 144 | + plt.close() |
| 145 | + |
| 146 | + |
| 147 | + # Print CO2 prices |
| 148 | + |
| 149 | + # for i, year in enumerate(years): |
| 150 | + # for scenario, decision in itertools.product(scenarios, decisions): |
| 151 | + |
| 152 | + # n = networks[year][scenario][decision] |
| 153 | + |
| 154 | + # print(f"CO2 price for {year}, {scenario}, {decision}: {n.global_constraints.loc["CO2Limit", "mu"] + n.global_constraints.loc["co2_limit-DE", "mu"]}") |
| 155 | + |
| 156 | + |
| 157 | + # Plot OPEX |
| 158 | + |
| 159 | + kwargs = { |
| 160 | + "groupby": ["bus", "carrier"], |
| 161 | + "at_port": True, |
| 162 | + "nice_names": False, |
| 163 | + } |
| 164 | + |
| 165 | + fig, axes = plt.subplots(nrows=3, ncols=1, figsize=(12, 18)) |
| 166 | + axes = axes.flatten() |
| 167 | + |
| 168 | + for i, year in enumerate(years): |
| 169 | + opex_comp = pd.DataFrame(columns=["_".join(tup) for tup in itertools.product(scenarios, decisions)]) |
| 170 | + |
| 171 | + # Collect OPEX for all scenario-decision combinations |
| 172 | + for scenario, decision in itertools.product(scenarios, decisions): |
| 173 | + n = networks[year][scenario][decision] |
| 174 | + |
| 175 | + opex = ( |
| 176 | + n.statistics.opex(**kwargs) |
| 177 | + .filter(like="DE") |
| 178 | + .groupby("carrier").sum() |
| 179 | + .multiply(1e-9) # to billion € |
| 180 | + ) |
| 181 | + opex_comp[f"{scenario}_{decision}"] = opex |
| 182 | + |
| 183 | + # Aggregate cost components with less than 0.1 (100 Mio €) as "Other" |
| 184 | + opex_comp_agg = aggregate_by_keywords(opex_comp, groups) |
| 185 | + small_rows = opex_comp_agg.abs().max(axis=1) < 0.1 |
| 186 | + other_row = opex_comp_agg[small_rows].sum(axis=0) |
| 187 | + opex_comp_agg = opex_comp_agg.loc[~small_rows] |
| 188 | + opex_comp_agg.loc['Other'] = other_row |
| 189 | + |
| 190 | + # Prepare labels with line breaks |
| 191 | + labels = [col.replace('_', '\n') for col in opex_comp_agg.columns] |
| 192 | + |
| 193 | + # Plot stacked bar |
| 194 | + ax = axes[i] |
| 195 | + bottom = np.zeros(len(opex_comp_agg.columns)) |
| 196 | + |
| 197 | + for tech in opex_comp_agg.index: |
| 198 | + values = opex_comp_agg.loc[tech].values |
| 199 | + ax.bar(labels, values, bottom=bottom, color=tech_colors.get(tech, '#333333'), label=tech) |
| 200 | + |
| 201 | + # Add numbers in the middle, except for 'Other' |
| 202 | + if tech != 'Other': |
| 203 | + for j, val in enumerate(values): |
| 204 | + if val > 0: # only if positive |
| 205 | + ax.text( |
| 206 | + j, |
| 207 | + bottom[j] + val/2, # middle of the segment |
| 208 | + f'{val:.2f}', |
| 209 | + ha='center', va='center', fontsize=8, color='white' |
| 210 | + ) |
| 211 | + |
| 212 | + bottom += values |
| 213 | + |
| 214 | + # Add total sum labels on top of bars |
| 215 | + totals = opex_comp_agg.sum(axis=0) |
| 216 | + for j, total in enumerate(totals): |
| 217 | + ax.text(j, total + total*0.02, f'{total:.2f}', ha='center', va='bottom', fontsize=10) |
| 218 | + |
| 219 | + # Adjust y-limit |
| 220 | + ax.set_ylim(0, max(totals)*1.08) |
| 221 | + ax.set_ylabel('OPEX [billion €]') |
| 222 | + ax.set_title(f'Stacked OPEX composition by technology, {year}') |
| 223 | + |
| 224 | + # Legend outside |
| 225 | + axes[-1].legend(loc='upper left', bbox_to_anchor=(1,1)) |
| 226 | + plt.savefig(snakemake.output[-1] + f"/opex_comp_de.png", bbox_inches="tight") |
| 227 | + plt.close() |
0 commit comments