Interconnector model

analyze_system.py

@@ -1,3 +1,4 @@
+from calendar import week
 from pathlib import Path
 
 import matplotlib.pyplot as plt
@@ -405,6 +406,75 @@ def plot_balance_week(n: pypsa.Network, output_dir: Path, season: str) -> None:
     plt.close(fig)
 
 
+
+def plot_balance_week2(n: pypsa.Network, output_dir: Path, season: str) -> None:
+    """
+    Plot the system balance for a representative week, with electricity demand
+    shown as a positive line.
+    """
+    season_weeks = get_representative_weeks(n.snapshots)
+
+    balance = n.statistics.energy_balance(aggregate_time=False)
+
+    # Aggregate by carrier
+    balance_by_carrier = balance.groupby(level="carrier").sum()
+
+    # Put time on x-axis
+    balance_by_carrier_t = balance_by_carrier.T
+
+    # Drop carriers that are all NaN
+    balance_by_carrier_t = balance_by_carrier_t.dropna(axis=1, how="all")
+
+    # Drop carriers that are zero everywhere
+    balance_by_carrier_t = balance_by_carrier_t.loc[
+        :, (balance_by_carrier_t.fillna(0).abs() > 0).any(axis=0)
+    ]
+
+    # Select representative week
+    week = balance_by_carrier_t.loc[season_weeks[season]].copy()
+
+    # Electricity demand as positive line
+    # Assumes electricity loads are attached to electricity buses
+    elec_loads = n.loads.index[n.loads.bus.map(n.buses.carrier) == "electricity"]
+    demand = n.loads_t.p_set[elec_loads].sum(axis=1).loc[season_weeks[season]]
+
+    relevant = ["gas", "offwind", "solar"]
+    week = week[relevant] if all(col in week.columns for col in relevant) else pd.DataFrame()
+
+    # Colors
+    colors = get_carrier_colors(n, week.columns)
+
+    # Styling
+    plt.rcParams.update({
+        "font.size": 13,
+        "axes.titlesize": 16,
+        "axes.labelsize": 16,
+        "legend.fontsize": 16,
+        "xtick.labelsize": 16,
+        "ytick.labelsize": 16,
+    })
+
+    fig, ax = plt.subplots(figsize=(13, 6))
+
+    # Positive and negative stacked areas
+    week.plot.area(ax=ax, stacked=True, color=colors)
+
+    # Demand line
+    demand.plot(ax=ax, color="black", linewidth=2, label="electricity demand")
+
+    ax.set_ylabel("Power [MW]")
+    ax.set_xlabel("")
+    ax.set_ylim(bottom=0)
+    ax.grid(axis="y", alpha=0.3)
+
+    # Put legend outside
+    ax.legend(loc="center left", bbox_to_anchor=(1.02, 0.5), frameon=False, title="")
+
+    fig.tight_layout()
+    fig.savefig(output_dir / f"balance_{season}_week.png", dpi=300, bbox_inches="tight")
+    plt.close(fig)
+
+
 def plot_duration_curves(n: pypsa.Network, output_dir: Path) -> None:
 
 
@@ -596,6 +666,50 @@ def plot_installed_capacity_by_generator(df: pd.DataFrame, output_path: str | No
     plt.show()
     plt.close(fig)
 
+def plot_capacity_factors_over_year(n: pypsa.Network, bus_carrier: str, output_dir: Path) -> None:
+    # Select generators at the given bus carrier
+    gens = n.generators.index[n.generators.bus.map(n.buses.carrier) == bus_carrier]
+
+    # Time series and capacities
+    p = n.generators_t.p[gens]
+    p_nom = n.generators.loc[gens, "p_nom_opt"]
+    weights = n.snapshot_weightings.generators
+
+    # Weighted generation
+    generation = p.multiply(weights, axis=0)
+
+    # 🔧 FIX: group by carrier (no axis=1)
+    generation = generation.T.groupby(n.generators.loc[gens, "carrier"]).sum().T
+
+    # Aggregate capacities by carrier
+    capacities = p_nom.groupby(n.generators.loc[gens, "carrier"]).sum()
+
+    # Aggregate over time (weekly or monthly)
+    generation = generation.resample("ME").sum()
+    hours = weights.resample("ME").sum()
+
+    # Capacity factor
+    cf = generation.divide(hours, axis=0)
+    cf = cf.divide(capacities, axis=1)
+
+    relevant = ["gas", "offwind", "solar"]
+    cf = cf[relevant]   
+
+    colors = get_carrier_colors(n, cf.columns)
+
+    # Plot
+    fig, ax = plt.subplots(figsize=(8, 4))
+    cf.plot(ax=ax, color=colors)
+
+    ax.set_ylabel("Capacity factor [-]")
+    ax.set_xlabel("")
+    ax.grid(axis="y", alpha=0.3)
+    ax.legend(loc="center left", bbox_to_anchor=(1, 0.5), title="")
+
+    fig.tight_layout()
+    fig.savefig(output_dir / "capacity_factors_over_year.png", dpi=300)
+    plt.close(fig)
+
 # Example usage:
 # plot_installed_capacity_by_generator(df, output_path="installed_capacity_by_generator.png")
 
@@ -661,13 +775,14 @@ def main() -> None:
     plot_annual_mix_from_balance(n, output_dir)
     plot_capacity_factors(n, bus_carrier, output_dir)
     plot_curtailment(n, bus_carrier, output_dir)
+    plot_capacity_factors_over_year(n, bus_carrier, output_dir)
 
 
     # Seasonal dispatch and duration curves
     plot_dispatch_week(dispatch_ts, season_weeks["winter"], "winter", output_dir)
     plot_dispatch_week(dispatch_ts, season_weeks["summer"], "summer", output_dir)
-    plot_balance_week(n, output_dir, "winter")
-    plot_balance_week(n, output_dir, "summer")
+    plot_balance_week2(n, output_dir, "winter")
+    plot_balance_week2(n, output_dir, "summer")
     plot_duration_curves(n, output_dir)
 
     # Interannual installed capacity by weather year

interconnectors_analyze.py

@@ -259,6 +259,189 @@ def plot_denmark_dispatch(n, folder):
     plt.savefig(folder / "denmark_dispatch_imports_battery_winter_week.png", dpi=300, bbox_inches='tight')
     plt.close()
 
+def plot_denmark_dispatch_strategy(n, folder):
+    """
+    Plot Denmark's dispatch strategy for the winter week with the highest
+    average Danish electricity demand.
+
+    The figure shows:
+    - Danish domestic generation by carrier (stacked area)
+    - imports from each neighbouring country (positive dashed lines)
+    - exports to each neighbouring country (negative dashed lines)
+    - Danish demand (black line)
+
+    Saves:
+        denmark_dispatch_strategy_winter_week.png
+    """
+
+    snapshots = n.snapshots
+
+    # Winter = December + January
+    winter_mask = (snapshots.month == 12) | (snapshots.month == 1)
+
+    # Denmark load
+    dk_load = n.loads_t.p_set["DK_electricity_demand"]
+    winter_load = dk_load[winter_mask]
+
+    # Find winter week with highest average load
+    weekly_avg_load = winter_load.resample("W").mean()
+    max_load_week_end = weekly_avg_load.idxmax()
+    week_start = max_load_week_end - pd.Timedelta(days=6)
+    week_end = max_load_week_end
+
+    week_index = dk_load.loc[week_start:week_end].index
+    week_load = dk_load.loc[week_index]
+
+    # -----------------------------
+    # 1. Denmark domestic generation by carrier
+    # -----------------------------
+    dk_generators = n.generators[n.generators.bus == "DK"]
+
+    generation_by_carrier = {}
+    for gen in dk_generators.index:
+        carrier = n.generators.at[gen, "carrier"]
+
+        if carrier not in generation_by_carrier:
+            generation_by_carrier[carrier] = pd.Series(0.0, index=week_index)
+
+        generation_by_carrier[carrier] = (
+            generation_by_carrier[carrier]
+            .add(n.generators_t.p[gen].loc[week_index], fill_value=0.0)
+        )
+
+    # Optional: include storage discharge/charge if your model uses links
+    # Battery charging is demand-like, battery discharging is supply-like
+    if "DK_battery_discharger" in n.links.index:
+        generation_by_carrier["battery_discharge"] = n.links_t.p1["DK_battery_discharger"].loc[week_index].clip(lower=0)
+
+    battery_charge = pd.Series(0.0, index=week_index)
+    if "DK_battery_charger" in n.links.index:
+        # charging consumes power from DK bus
+        # depending on model convention p0 may already be positive when consuming
+        battery_charge = n.links_t.p0["DK_battery_charger"].loc[week_index].clip(lower=0)
+
+    # -----------------------------
+    # 2. Imports/exports by neighbour
+    # -----------------------------
+    dk_lines = n.lines[(n.lines.bus0 == "DK") | (n.lines.bus1 == "DK")]
+
+    exchanges = {}
+
+    for line in dk_lines.index:
+        bus0 = n.lines.at[line, "bus0"]
+        bus1 = n.lines.at[line, "bus1"]
+
+        neighbour = bus1 if bus0 == "DK" else bus0
+        flow = n.lines_t.p0[line].loc[week_index]
+
+        # PyPSA convention:
+        # p0 > 0 means power flows from bus0 -> bus1
+        # Convert to "positive = import to DK, negative = export from DK"
+        if bus0 == "DK":
+            dk_exchange = -flow
+        else:
+            dk_exchange = flow
+
+        if neighbour not in exchanges:
+            exchanges[neighbour] = pd.Series(0.0, index=week_index)
+
+        exchanges[neighbour] = exchanges[neighbour].add(dk_exchange, fill_value=0.0)
+
+    # -----------------------------
+    # 3. Plot (two panels)
+    # -----------------------------
+
+    preferred_order = [
+    "solar", "onwind", "offwind",
+    "gas", "coal", "nuclear",
+    "battery_discharge"
+    ]
+
+    remaining = [c for c in generation_by_carrier if c not in preferred_order]
+
+    carrier_order = [c for c in preferred_order if c in generation_by_carrier] + remaining
+
+    fig, (ax1, ax2) = plt.subplots(
+        2, 1,
+        figsize=(14, 9),
+        sharex=True,
+        gridspec_kw={"height_ratios": [3, 1]}
+    )
+
+    # -----------------------------
+    # TOP PANEL: DK generation
+    # -----------------------------
+    colors = {
+        "solar": "#ffd92f",
+        "onwind": "#1b9e77",
+        "offwind": "#377eb8",
+        "gas": "#e41a1c",
+        "coal": "#4d4d4d",
+        "nuclear": "#984ea3",
+        "battery_discharge": "#ff7f00",
+    }
+
+    bottom = np.zeros(len(week_index))
+
+    for carrier in carrier_order:
+        series = generation_by_carrier[carrier].fillna(0).values
+
+        ax1.fill_between(
+            week_index,
+            bottom,
+            bottom + series,
+            label=f"DK {carrier}",
+            color=colors.get(carrier, None),
+            alpha=0.8
+        )
+
+        bottom += series
+
+    # Demand
+    ax1.plot(
+        week_index,
+        week_load,
+        color="black",
+        linewidth=2.3,
+        label="DK demand"
+    )
+
+    ax1.set_ylabel("Power [MW]")
+    ax1.set_title("Denmark dispatch strategy during peak winter week")
+    ax1.legend(loc="upper left", bbox_to_anchor=(1.01, 1))
+    ax1.grid(True, alpha=0.3)
+
+    # -----------------------------
+    # BOTTOM PANEL: Exchanges
+    # -----------------------------
+    exchange_colors = {
+        "DE": "#2ca02c",
+        "SE": "#ff7f0e",
+        "NO": "#1f77b4"
+    }
+
+    for neighbour, series in exchanges.items():
+        ax2.plot(
+            week_index,
+            series,
+            linewidth=2.5,
+            label=f"{neighbour}",
+            color=exchange_colors.get(neighbour, None)
+        )
+
+    ax2.axhline(0, color="black", linewidth=1)
+    ax2.set_ylabel("Import / Export [MW]")
+    ax2.set_xlabel("Time")
+    ax2.legend(title="Exchange", loc="upper left", bbox_to_anchor=(1.01, 1))
+    ax2.grid(True, alpha=0.3)
+
+    plt.xticks(rotation=45)
+    plt.tight_layout()
+
+    outfile = folder / "denmark_dispatch_strategy_winter_week.png"
+    plt.savefig(outfile, dpi=300, bbox_inches="tight")
+    plt.close()
+
 
 def main():
     """Main function to run the analysis."""
@@ -282,6 +465,7 @@ def main():
     plot_topology(n, folder)
     plot_line_loading_duration_curves(n, folder)
     plot_denmark_dispatch(n, folder)
+    plot_denmark_dispatch_strategy(n, folder)
 
     print(f"Analysis complete. Files saved in {folder}")