3-Node Capacity Expansion¶

This example builds on the single-node capacity expansion example and extends it to three regions in Australia (New South Wales, Victoria, and South Australia), which can be connected via transmission links. The purpose is to illustrate how to set up a capacity expansion model with multiple interconnected regions, where the model can optimize for trade-offs between generation capacity, storage, and transmission infrastructure.

In this first block, we import the necessary libraries and set up some basic functions (e.g. annuity factors) and constants (e.g. temporal resolution and solver):

In [2]:

import cartopy.crs as ccrs
import geopandas as gpd
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

import pypsa
from pypsa.costs import annuity

RESOLUTION = 3  # 3-hourly

import cartopy.crs as ccrs import geopandas as gpd import matplotlib.pyplot as plt import numpy as np import pandas as pd import pypsa from pypsa.costs import annuity RESOLUTION = 3 # 3-hourly

In a fresh PyPSA network instance, we set the coordinates of each state's capital (Sydney, Melbourne, and Adelaide) as the coordinates of the electricity buses.

In [3]:

n = pypsa.Network()

REGIONS = pd.Index(["VIC", "SA", "NSW"])
LAT = [-37.81, -34.93, -33.87]
LON = [145.01, 138.63, 151.20]
n.add("Bus", REGIONS, x=LON, y=LAT);

n = pypsa.Network() REGIONS = pd.Index(["VIC", "SA", "NSW"]) LAT = [-37.81, -34.93, -33.87] LON = [145.01, 138.63, 151.20] n.add("Bus", REGIONS, x=LON, y=LAT);

We also define a range of carriers (i.e. energy carriers or technologies) that we will use in the model. This will be useful later when plotting optimisation results.

In [4]:

CARRIERS = {
    "solar": "gold",
    "wind": "steelblue",
    "load shedding": "indianred",
    "battery charger": "lightgreen",
    "battery discharger": "lightgreen",
    "battery storage": "grey",
    "battery": "grey",
    "electrolysis": "violet",
    "turbine": "violet",
    "hydrogen storage": "orchid",
    "hydrogen": "orchid",
    "AC": "black",
    "HVDC": "lightseagreen",
    "load": "slategrey",
}

n.add(
    "Carrier",
    CARRIERS.keys(),
    color=CARRIERS.values(),
);

CARRIERS = { "solar": "gold", "wind": "steelblue", "load shedding": "indianred", "battery charger": "lightgreen", "battery discharger": "lightgreen", "battery storage": "grey", "battery": "grey", "electrolysis": "violet", "turbine": "violet", "hydrogen storage": "orchid", "hydrogen": "orchid", "AC": "black", "HVDC": "lightseagreen", "load": "slategrey", } n.add( "Carrier", CARRIERS.keys(), color=CARRIERS.values(), );

As data, we have given:

• australia-example-p_set.csv - a time series of electricity demand in each region in 2011 from the Australian Energy Market Operator (AEMO)
• australia-example-p_max_pu.csv - a time series of capacity factors in 2011 for onshore wind and solar PV for each region from model.energy
• australia-example-shapes.geojson - a geographical representation of the three regions from the GeoBoundaries project

The timestamps in the datasets are given in UTC.

The datasets have been obtained from the following code:

import pandas as pd
import geopandas as gpd

REGIONS = ["VIC", "SA", "NSW"]

# 1. Capacity factors
capacity_factors = {
    "VIC": "https://model.energy/data/time-series-c3ec7795e4d52617defed46495899b50.csv",
    "SA": "https://model.energy/data/time-series-efcb605660bd6f9db75fc7c5b3fd5379.csv",
    "NSW": "https://model.energy/data/time-series-ef16f3da1c24442d7efac9d7c4a0ce05.csv",
}

p_max_pu = pd.concat({
    c: pd.read_csv(
        url,
        index_col=0,
        parse_dates=True,
    )
    for c, url in capacity_factors.items()
}, axis=1)

p_max_pu.to_csv("australia-example-p_max_pu.csv")

# 2. Electricity demand
so = {'User-Agent': 'Mozilla/5.0'}
year_months = [f"2011{m:02d}" for m in range(1, 13)] + ["201201"]  # All months in 2011 + Jan 2012

p_set = pd.concat([
    pd.concat([
        pd.read_csv(
            f"https://aemo.com.au/aemo/data/nem/priceanddemand/PRICE_AND_DEMAND_{ym}_{region}1.csv",
            storage_options=so,
            index_col=1,
            parse_dates=True
        )["TOTALDEMAND"].rename(region)
        for ym in year_months
    ])
    for region in REGIONS
], axis=1)

p_set = p_set.resample("1h").mean().ffill(inplace=True)

p_set = p_set.tz_localize(
    "Australia/NSW", ambiguous="NaT", nonexistent="shift_backward"
).tz_convert("UTC").tz_localize(None).reindex(index=p_max_pu.index)

p_set.to_csv("australia-example-p_set.csv")

# 3. Geographical shapes
url = "https://media.githubusercontent.com/media/wmgeolab/geoBoundaries/bdfb316b1fdcac1473051979152cac0943c549fa/releaseData/gbOpen/AUS/ADM1/geoBoundaries-AUS-ADM1-all.zip"

gdf = gpd.read_file(url, layer="geoBoundaries-AUS-ADM1_simplified")
gdf.index = gdf.shapeISO.str.split("-").str[1]
gdf = gdf.loc[REGIONS]
gdf.to_file("australia-example-shapes.geojson", driver="GeoJSON")

We load this data and set the snapshots of the network to the timestamps of the time series data.

In [5]:

url = "https://tubcloud.tu-berlin.de/s/oPQbAebrciFBZP2/download/australia-example-p_max_pu.csv"
p_max_pu = pd.read_csv(url, index_col=0, parse_dates=True, header=[0, 1])
display(p_max_pu.head(3))

url = "https://tubcloud.tu-berlin.de/s/8C6d7z3HxE7yZi9/download/australia-example-p_set.csv"
p_set = pd.read_csv(url, index_col=0, parse_dates=True)
display(p_set.head(3))

n.set_snapshots(p_max_pu.index)

url = "https://tubcloud.tu-berlin.de/s/oPQbAebrciFBZP2/download/australia-example-p_max_pu.csv" p_max_pu = pd.read_csv(url, index_col=0, parse_dates=True, header=[0, 1]) display(p_max_pu.head(3)) url = "https://tubcloud.tu-berlin.de/s/8C6d7z3HxE7yZi9/download/australia-example-p_set.csv" p_set = pd.read_csv(url, index_col=0, parse_dates=True) display(p_set.head(3)) n.set_snapshots(p_max_pu.index)

	VIC		SA		NSW
	onwind	solar	onwind	solar	onwind	solar
2011-01-01 00:00:00	0.214	0.504	0.239	0.401	0.141	0.515
2011-01-01 01:00:00	0.197	0.588	0.244	0.504	0.098	0.597
2011-01-01 02:00:00	0.215	0.613	0.232	0.568	0.068	0.644

	VIC	SA	NSW
snapshot
2011-01-01 00:00:00	5113.930	1358.385	9357.435
2011-01-01 01:00:00	5147.425	1383.295	9635.540
2011-01-01 02:00:00	5079.465	1387.215	9832.175

Then, we add the loads for each region:

In [6]:

n.add("Load", REGIONS, suffix=" load", bus=REGIONS, p_set=p_set, carrier="load");

n.add("Load", REGIONS, suffix=" load", bus=REGIONS, p_set=p_set, carrier="load");

For the wind and solar generators, we pass the capacity factors as p_max_pu and add the annuitized capital costs. Here, we assume a capital cost of 2000 AUD/kW for wind and 700 AUD/kW for solar, with a discount rate of 5% and a lifetime of 20 years. Note again that the capacity is given in MW, so we multiply the capital costs by 1000 to convert them to kW.

As a trick for modelling curtailment, we add a small marginal_cost to the wind generators to nudge those to be curtailed first before solar.

In [7]:

p_max_pu_wind = p_max_pu.xs("onwind", level=1, axis=1).rename(
    columns=lambda s: s + " wind"
)

n.add(
    "Generator",
    p_max_pu_wind.columns,
    bus=REGIONS,
    p_max_pu=p_max_pu_wind,
    p_nom_extendable=True,
    capital_cost=annuity(0.05, 30) * 2_000_000,
    marginal_cost=0.5,
    carrier="wind",
)

p_max_pu_solar = p_max_pu.xs("solar", level=1, axis=1).rename(
    columns=lambda s: s + " solar"
)

n.add(
    "Generator",
    p_max_pu_solar.columns,
    bus=REGIONS,
    p_max_pu=p_max_pu_solar,
    p_nom_extendable=True,
    capital_cost=annuity(0.05, 30) * 700_000,
    carrier="solar",
);

p_max_pu_wind = p_max_pu.xs("onwind", level=1, axis=1).rename( columns=lambda s: s + " wind" ) n.add( "Generator", p_max_pu_wind.columns, bus=REGIONS, p_max_pu=p_max_pu_wind, p_nom_extendable=True, capital_cost=annuity(0.05, 30) * 2_000_000, marginal_cost=0.5, carrier="wind", ) p_max_pu_solar = p_max_pu.xs("solar", level=1, axis=1).rename( columns=lambda s: s + " solar" ) n.add( "Generator", p_max_pu_solar.columns, bus=REGIONS, p_max_pu=p_max_pu_solar, p_nom_extendable=True, capital_cost=annuity(0.05, 30) * 700_000, carrier="solar", );

We also add a generator with a very high marginal cost to represent load shedding with a value of lost load (VoLL) of 3000 AUD/MWh.

In [8]:

n.add(
    "Generator",
    REGIONS,
    suffix=" load shedding",
    bus=REGIONS,
    p_nom=p_set.max(),
    marginal_cost=3000,
    carrier="load shedding",
);

n.add( "Generator", REGIONS, suffix=" load shedding", bus=REGIONS, p_nom=p_set.max(), marginal_cost=3000, carrier="load shedding", );

To save some computation time, we only sample every third snapshot, which corresponds to a temporal resolution of 3 hours. Note that the snapshot weightings have to be adjusted accordingly.

Let's give this basic model a first spin:

In [9]:

n.set_snapshots(n.snapshots[::RESOLUTION])
n.snapshot_weightings.loc[:, :] = RESOLUTION

n.optimize(solver_name="gurobi")

n.set_snapshots(n.snapshots[::RESOLUTION]) n.snapshot_weightings.loc[:, :] = RESOLUTION n.optimize(solver_name="gurobi")

/tmp/ipykernel_3237/3103030024.py:4: FutureWarning: The default value of `include_objective_constant` will change from True to False in version 2.0. Set `include_objective_constant` explicitly to suppress this warning. Using False improves LP numerical conditioning by not including the objective constant as a variable.
  n.optimize(solver_name="gurobi")

INFO:linopy.model: Solve problem using Gurobi solver

INFO:linopy.model:Solver options:
 - log_to_console: False

INFO:linopy.io:Writing objective.

Writing constraints.:   0%|          | 0/6 [00:00<?, ?it/s]

Writing constraints.: 100%|██████████| 6/6 [00:00<00:00, 110.77it/s]

Writing continuous variables.:   0%|          | 0/2 [00:00<?, ?it/s]

Writing continuous variables.: 100%|██████████| 2/2 [00:00<00:00, 254.74it/s]

INFO:linopy.io: Writing time: 0.08s

Set parameter WLSAccessID

Set parameter WLSSecret

Set parameter LicenseID to value 2537914

Academic license 2537914 - for non-commercial use only - registered to l.___@tu-berlin.de

Read LP format model from file /tmp/linopy-problem-j6g3414r.lp

Reading time = 0.07 seconds

obj: 61326 rows, 26286 columns, 92010 nonzeros

Set parameter LogToConsole to value 0

Warning: environment still referenced so free is deferred (Continue to use WLS)

INFO:linopy.constants: Optimization successful: 
Status: ok
Termination condition: optimal
Solution: 26286 primals, 61326 duals
Objective: 3.16e+10
Solver model: available
Solver message: 2

INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-fix-p-lower, Generator-fix-p-upper, Generator-ext-p-lower, Generator-ext-p-upper were not assigned to the network.

Out[9]:

('ok', 'optimal')

In this first run, every region is isolated (no transmission) and has no storage (no batteries or hydrogen storage).

Around 2.5% of the load is shed; roughly 50 TWh come from solar and 80 TWh from wind. The average cost of 226.57 AUD/MWh is quite high.

In [10]:

display(n.statistics.energy_balance().div(1e6).round(2).sort_values())  # TWh

average_cost = (
    (n.statistics.capex().sum() + n.statistics.opex().sum())
    / n.loads_t.p_set.sum().sum()
    / RESOLUTION
)

display(f"Average cost: {average_cost:.2f} AUD/MWh")

display(n.statistics.energy_balance().div(1e6).round(2).sort_values()) # TWh average_cost = ( (n.statistics.capex().sum() + n.statistics.opex().sum()) / n.loads_t.p_set.sum().sum() / RESOLUTION ) display(f"Average cost: {average_cost:.2f} AUD/MWh")

component  carrier        bus_carrier
Load       load           AC            -139.62
Generator  load shedding  AC               3.25
           solar          AC              52.65
           wind           AC              83.73
dtype: float64

'Average cost: 226.57 AUD/MWh'

The statistics and plotting functionality of PyPSA can be used to create a more detailed overview of the results on a map. For example, we can plot the energy balance and energy mix of each region as pie charts, where the upper half shows the generation and the lower half the consumption. We can also add other geographical layers from other objects to the map, such as the geographical shapes of the regions colored by the average locational marginal price (LMP).

In [11]:

# Create an empty figure with a Cartopy projection
crs = ccrs.PlateCarree()
fig, ax = plt.subplots(figsize=(10, 6), subplot_kw={"projection": crs})

# Use the energy balance statistics to prepare the bus sizes and plot the network
bus_size = n.statistics.energy_balance(groupby=["bus", "carrier"]).droplevel(
    "component"
)
n.plot(ax=ax, bus_size=bus_size / 4e7, margin=0.75, bus_split_circle=True)

# Load the shapefiles and add them to the plot colored by average LMP
prices = n.buses_t.marginal_price.mean()
gdf = gpd.read_file(
    "https://tubcloud.tu-berlin.de/s/4n9PegitETESBGR/download/australia-example-shapes.geojson"
).set_index("index")
gdf.to_crs(crs).plot(ax=ax, column=prices, cmap="RdYlGn_r", vmin=90, vmax=240)

# Add a legend for the LMP color scale
norm = plt.Normalize(vmin=90, vmax=240)
sm = plt.cm.ScalarMappable(cmap="RdYlGn_r", norm=norm)
fig.colorbar(sm, ax=ax, label="LMP [AUD/MWh]", shrink=0.4)

# Create an empty figure with a Cartopy projection crs = ccrs.PlateCarree() fig, ax = plt.subplots(figsize=(10, 6), subplot_kw={"projection": crs}) # Use the energy balance statistics to prepare the bus sizes and plot the network bus_size = n.statistics.energy_balance(groupby=["bus", "carrier"]).droplevel( "component" ) n.plot(ax=ax, bus_size=bus_size / 4e7, margin=0.75, bus_split_circle=True) # Load the shapefiles and add them to the plot colored by average LMP prices = n.buses_t.marginal_price.mean() gdf = gpd.read_file( "https://tubcloud.tu-berlin.de/s/4n9PegitETESBGR/download/australia-example-shapes.geojson" ).set_index("index") gdf.to_crs(crs).plot(ax=ax, column=prices, cmap="RdYlGn_r", vmin=90, vmax=240) # Add a legend for the LMP color scale norm = plt.Normalize(vmin=90, vmax=240) sm = plt.cm.ScalarMappable(cmap="RdYlGn_r", norm=norm) fig.colorbar(sm, ax=ax, label="LMP [AUD/MWh]", shrink=0.4)

/home/docs/checkouts/readthedocs.org/user_builds/pypsa/envs/latest/lib/python3.13/site-packages/cartopy/io/__init__.py:242: DownloadWarning: Downloading: https://naturalearth.s3.amazonaws.com/50m_cultural/ne_50m_admin_0_boundary_lines_land.zip
  warnings.warn(f'Downloading: {url}', DownloadWarning)

/home/docs/checkouts/readthedocs.org/user_builds/pypsa/envs/latest/lib/python3.13/site-packages/cartopy/io/__init__.py:242: DownloadWarning: Downloading: https://naturalearth.s3.amazonaws.com/50m_physical/ne_50m_coastline.zip
  warnings.warn(f'Downloading: {url}', DownloadWarning)

Out[11]:

<matplotlib.colorbar.Colorbar at 0x7a0b22c2b0e0>

Now, we add some storage options to each model region:

• For battery storage, we separately model the battery inverter (charger and discharger) and the battery storage itself. This has the advantage that the power-to-energy ratio can be optimised.
• For hydrogen storage, we model the electrolyser, the hydrogen storage, and the turbine for re-electrification separately, so that these components can be independently sized.

In [12]:

n.add("Bus", REGIONS, suffix=" hydrogen", carrier="hydrogen", x=LON, y=LAT)

n.add(
    "Link",
    REGIONS,
    suffix=" electrolysis",
    bus0=REGIONS,
    bus1=pd.Index(REGIONS) + " hydrogen",
    carrier="electrolysis",
    p_nom_extendable=True,
    efficiency=0.7,
    capital_cost=annuity(0.05, 25) * 2_500_000,
)

n.add(
    "Link",
    REGIONS,
    suffix=" turbine",
    bus0=pd.Index(REGIONS) + " hydrogen",
    bus1=REGIONS,
    carrier="turbine",
    p_nom_extendable=True,
    efficiency=0.4,
    capital_cost=annuity(0.05, 25) * 2_000_000,
)

n.add(
    "Store",
    REGIONS,
    suffix=" hydrogen storage",
    bus=pd.Index(REGIONS) + " hydrogen",
    carrier="hydrogen storage",
    capital_cost=annuity(0.05, 30) * 80,
    e_nom_extendable=True,
    e_cyclic=True,
)

n.add("Bus", REGIONS, suffix=" battery", carrier="battery", x=LON, y=LAT)

n.add(
    "Link",
    REGIONS,
    suffix=" battery charger",
    bus0=REGIONS,
    bus1=REGIONS + " battery",
    carrier="battery charger",
    p_nom_extendable=True,
    efficiency=0.95,
    capital_cost=annuity(0.05, 10) * 300_000,
)

n.add(
    "Link",
    REGIONS,
    suffix=" battery discharger",
    bus0=REGIONS + " battery",
    bus1=REGIONS,
    carrier="battery discharger",
    p_nom_extendable=True,
    efficiency=0.95,
)

n.add(
    "Store",
    REGIONS,
    suffix=" battery storage",
    bus=REGIONS + " battery",
    carrier="battery storage",
    capital_cost=annuity(0.05, 25) * 250_000,
    e_nom_extendable=True,
    e_cyclic=True,
);

n.add("Bus", REGIONS, suffix=" hydrogen", carrier="hydrogen", x=LON, y=LAT) n.add( "Link", REGIONS, suffix=" electrolysis", bus0=REGIONS, bus1=pd.Index(REGIONS) + " hydrogen", carrier="electrolysis", p_nom_extendable=True, efficiency=0.7, capital_cost=annuity(0.05, 25) * 2_500_000, ) n.add( "Link", REGIONS, suffix=" turbine", bus0=pd.Index(REGIONS) + " hydrogen", bus1=REGIONS, carrier="turbine", p_nom_extendable=True, efficiency=0.4, capital_cost=annuity(0.05, 25) * 2_000_000, ) n.add( "Store", REGIONS, suffix=" hydrogen storage", bus=pd.Index(REGIONS) + " hydrogen", carrier="hydrogen storage", capital_cost=annuity(0.05, 30) * 80, e_nom_extendable=True, e_cyclic=True, ) n.add("Bus", REGIONS, suffix=" battery", carrier="battery", x=LON, y=LAT) n.add( "Link", REGIONS, suffix=" battery charger", bus0=REGIONS, bus1=REGIONS + " battery", carrier="battery charger", p_nom_extendable=True, efficiency=0.95, capital_cost=annuity(0.05, 10) * 300_000, ) n.add( "Link", REGIONS, suffix=" battery discharger", bus0=REGIONS + " battery", bus1=REGIONS, carrier="battery discharger", p_nom_extendable=True, efficiency=0.95, ) n.add( "Store", REGIONS, suffix=" battery storage", bus=REGIONS + " battery", carrier="battery storage", capital_cost=annuity(0.05, 25) * 250_000, e_nom_extendable=True, e_cyclic=True, );

For the battery storage, we need to take special care of the inverter. As the same component can be used for charging and discharging, the capacities of our battery charger and discharger component must be linked. This is not done automatically by PyPSA, so we have to add an extra constraint, which we can pass as extra_functionality to the optimisation.

In [13]:

def battery_constraint(n: pypsa.Network, sns: pd.Index) -> None:
    """Constraint to ensure that the nominal capacity of battery chargers and dischargers are in a fixed ratio."""
    dischargers_i = n.links[n.links.index.str.contains(" discharger")].index
    chargers_i = n.links[n.links.index.str.contains(" charger")].index

    eff = n.links.efficiency[dischargers_i].values
    lhs = n.model["Link-p_nom"].loc[chargers_i]
    rhs = n.model["Link-p_nom"].loc[dischargers_i] * eff

    n.model.add_constraints(lhs == rhs, name="Link-charger_ratio")


n.optimize(extra_functionality=battery_constraint, solver_name="gurobi")

def battery_constraint(n: pypsa.Network, sns: pd.Index) -> None: """Constraint to ensure that the nominal capacity of battery chargers and dischargers are in a fixed ratio.""" dischargers_i = n.links[n.links.index.str.contains(" discharger")].index chargers_i = n.links[n.links.index.str.contains(" charger")].index eff = n.links.efficiency[dischargers_i].values lhs = n.model["Link-p_nom"].loc[chargers_i] rhs = n.model["Link-p_nom"].loc[dischargers_i] * eff n.model.add_constraints(lhs == rhs, name="Link-charger_ratio") n.optimize(extra_functionality=battery_constraint, solver_name="gurobi")

/tmp/ipykernel_3237/2124492770.py:13: FutureWarning: The default value of `include_objective_constant` will change from True to False in version 2.0. Set `include_objective_constant` explicitly to suppress this warning. Using False improves LP numerical conditioning by not including the objective constant as a variable.
  n.optimize(extra_functionality=battery_constraint, solver_name="gurobi")

INFO:linopy.model: Solve problem using Gurobi solver

INFO:linopy.model:Solver options:
 - log_to_console: False

INFO:linopy.io:Writing objective.

Writing constraints.:   0%|          | 0/14 [00:00<?, ?it/s]

Writing constraints.:  79%|███████▊  | 11/14 [00:00<00:00, 97.69it/s]

Writing constraints.: 100%|██████████| 14/14 [00:00<00:00, 79.97it/s]

Writing continuous variables.:   0%|          | 0/7 [00:00<?, ?it/s]

Writing continuous variables.: 100%|██████████| 7/7 [00:00<00:00, 259.42it/s]

INFO:linopy.io: Writing time: 0.21s

Set parameter WLSAccessID

Set parameter WLSSecret

Set parameter LicenseID to value 2537914

Academic license 2537914 - for non-commercial use only - registered to l.___@tu-berlin.de

Read LP format model from file /tmp/linopy-problem-_az7x9iv.lp

Reading time = 0.23 seconds

obj: 201507 rows, 96384 columns, 389874 nonzeros

Set parameter LogToConsole to value 0

Warning: environment still referenced so free is deferred (Continue to use WLS)

INFO:linopy.constants: Optimization successful: 
Status: ok
Termination condition: optimal
Solution: 96384 primals, 201507 duals
Objective: 1.40e+10
Solver model: available
Solver message: 2

INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-fix-p-lower, Generator-fix-p-upper, Generator-ext-p-lower, Generator-ext-p-upper, Link-ext-p-lower, Link-ext-p-upper, Store-ext-e-lower, Store-ext-e-upper, Store-energy_balance were not assigned to the network.

Out[13]:

('ok', 'optimal')

A look at the new energy balance shows use of both battery and hydrogen storage, as well as a shift from wind to solar generation compared to the previous run. In fact, overall, we see more generation from wind and solar than the load, in order to compansate for the losses in the storage systems. There is a difference of around 5 TWh between power consumption and injection of the battery and a 15 TWh difference between electrolyser electricity consumption and turbine electricity production. Yet, overall, average costs are reduced drastically to 100.51 AUD/MWh.

In [14]:

display(
    n.statistics.energy_balance(bus_carrier="AC").div(1e6).round(2).sort_values()
)  # TWh

average_cost = (
    (n.statistics.capex().sum() + n.statistics.opex().sum())
    / n.loads_t.p_set.sum().sum()
    / RESOLUTION
)

print(f"Average cost: {average_cost:.2f} AUD/MWh")

display( n.statistics.energy_balance(bus_carrier="AC").div(1e6).round(2).sort_values() ) # TWh average_cost = ( (n.statistics.capex().sum() + n.statistics.opex().sum()) / n.loads_t.p_set.sum().sum() / RESOLUTION ) print(f"Average cost: {average_cost:.2f} AUD/MWh")

component  carrier             bus_carrier
Load       load                AC            -139.62
Link       battery charger     AC             -51.11
           electrolysis        AC             -21.64
Generator  load shedding       AC               0.09
Link       turbine             AC               6.06
Generator  wind                AC              45.80
Link       battery discharger  AC              46.12
Generator  solar               AC             114.29
dtype: float64

Average cost: 100.51 AUD/MWh

Again, we can show these results spatially disagreggated on a map. The upper and lower half circles are still equal because no imbalances are possible due to the lack of transmission links. However, we can now see green elements for battery and purple elements for hydrogen storage, in addition to wind and solar generation.

In [15]:

bus_size = (
    n.statistics.energy_balance(groupby=["bus", "carrier"])
    .droplevel("component")
    .loc[REGIONS]
)

n.plot(bus_size=bus_size / 4e7, margin=0.75, bus_split_circle=True);

bus_size = ( n.statistics.energy_balance(groupby=["bus", "carrier"]) .droplevel("component") .loc[REGIONS] ) n.plot(bus_size=bus_size / 4e7, margin=0.75, bus_split_circle=True);

Of course, there are other plots we can create from the optimised model, for instance, a price duration curve for Victoria. For this we sort the time series of bus marginal prices in descending order and reindex to a relative axis from 0% to 100%:

In [16]:

pdc = (
    n.buses_t.marginal_price["VIC"].sort_values(ascending=False).reset_index(drop=True)
)
pdc.index = np.arange(0, 100, 100 / len(pdc.index))
pdc.plot(
    ylim=[0, 1000], xlim=[0, 100], xlabel="Share of time [%]", ylabel="Price [AUD/MWh]"
)

pdc = ( n.buses_t.marginal_price["VIC"].sort_values(ascending=False).reset_index(drop=True) ) pdc.index = np.arange(0, 100, 100 / len(pdc.index)) pdc.plot( ylim=[0, 1000], xlim=[0, 100], xlabel="Share of time [%]", ylabel="Price [AUD/MWh]" )

Out[16]:

<Axes: xlabel='Share of time [%]', ylabel='Price [AUD/MWh]'>

As a final modelling step, we will now add transmission links between the regions. By modelling the interconnectors as bidirectional links, we will assume that they are controllable point-to-point HVDC connections (rather than lines subject to Kirchhoff's Voltage Law). They are also assumed to be lossless, otherwise a similar workaround with two unidirectional links as for battery inverters has to be implemented. The cost of the transmission links is length-dependent and set to 1000 AUD/MW/km.

In [17]:

connections = [
    ("VIC", "NSW", 700),  # km
    ("VIC", "SA", 650),  # km
    ("NSW", "SA", 1150),  # km
]

for bus0, bus1, length in connections:
    n.add(
        "Link",
        f"{bus0}-{bus1}",
        bus0=bus0,
        bus1=bus1,
        carrier="HVDC",
        p_nom_extendable=True,
        capital_cost=annuity(0.05, 40) * 1_000 * length,
        p_min_pu=-1,  # bidirectional
    )

connections = [ ("VIC", "NSW", 700), # km ("VIC", "SA", 650), # km ("NSW", "SA", 1150), # km ] for bus0, bus1, length in connections: n.add( "Link", f"{bus0}-{bus1}", bus0=bus0, bus1=bus1, carrier="HVDC", p_nom_extendable=True, capital_cost=annuity(0.05, 40) * 1_000 * length, p_min_pu=-1, # bidirectional )

In [18]:

n.optimize(extra_functionality=battery_constraint, solver_name="gurobi")

n.optimize(extra_functionality=battery_constraint, solver_name="gurobi")

/tmp/ipykernel_3237/1825293787.py:1: FutureWarning: The default value of `include_objective_constant` will change from True to False in version 2.0. Set `include_objective_constant` explicitly to suppress this warning. Using False improves LP numerical conditioning by not including the objective constant as a variable.
  n.optimize(extra_functionality=battery_constraint, solver_name="gurobi")

INFO:linopy.model: Solve problem using Gurobi solver

INFO:linopy.model:Solver options:
 - log_to_console: False

INFO:linopy.io:Writing objective.

Writing constraints.:   0%|          | 0/14 [00:00<?, ?it/s]

Writing constraints.:  71%|███████▏  | 10/14 [00:00<00:00, 96.08it/s]

Writing constraints.: 100%|██████████| 14/14 [00:00<00:00, 73.40it/s]

Writing continuous variables.:   0%|          | 0/7 [00:00<?, ?it/s]

Writing continuous variables.: 100%|██████████| 7/7 [00:00<00:00, 238.42it/s]

INFO:linopy.io: Writing time: 0.23s

Set parameter WLSAccessID

Set parameter WLSSecret

Set parameter LicenseID to value 2537914

Academic license 2537914 - for non-commercial use only - registered to l.___@tu-berlin.de

Read LP format model from file /tmp/linopy-problem-23s36tmi.lp

Reading time = 0.25 seconds

obj: 219030 rows, 105147 columns, 442437 nonzeros

Set parameter LogToConsole to value 0

Warning: environment still referenced so free is deferred (Continue to use WLS)

INFO:linopy.constants: Optimization successful: 
Status: ok
Termination condition: optimal
Solution: 105147 primals, 219030 duals
Objective: 1.30e+10
Solver model: available
Solver message: 2

INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-fix-p-lower, Generator-fix-p-upper, Generator-ext-p-lower, Generator-ext-p-upper, Link-ext-p-lower, Link-ext-p-upper, Store-ext-e-lower, Store-ext-e-upper, Store-energy_balance were not assigned to the network.

Out[18]:

('ok', 'optimal')

With transmission enabled, we can see that the average costs are reduced further to 93.00 AUD/MWh. This is achieved as the interconnection facilitates the integration of wind variability across regions, reducing curtailment and the some need for hydrogen storage. The energy mix also shifts slightly from solar towards wind.

In [19]:

display(
    n.statistics.energy_balance(bus_carrier="AC").div(1e6).round(2).sort_values()
)  # TWh

average_cost = (
    (n.statistics.capex().sum() + n.statistics.opex().sum())
    / n.loads_t.p_set.sum().sum()
    / RESOLUTION
)

display(f"Average cost: {average_cost:.2f} AUD/MWh")

display( n.statistics.energy_balance(bus_carrier="AC").div(1e6).round(2).sort_values() ) # TWh average_cost = ( (n.statistics.capex().sum() + n.statistics.opex().sum()) / n.loads_t.p_set.sum().sum() / RESOLUTION ) display(f"Average cost: {average_cost:.2f} AUD/MWh")

component  carrier             bus_carrier
Load       load                AC            -139.62
Link       battery charger     AC             -51.55
           electrolysis        AC             -18.01
Generator  load shedding       AC               0.05
Link       turbine             AC               5.04
           battery discharger  AC              46.52
Generator  wind                AC              46.55
           solar               AC             111.02
dtype: float64

'Average cost: 93.00 AUD/MWh'

Now, we can add further elements to the map, such as the interconnector capacities, the net flow direction as well as the average loading of the interconnectors, which show how Victoria becomes a net importer of electricity from New South Wales and South Australia.

In [20]:

bus_size = (
    n.statistics.energy_balance(bus_carrier="AC", groupby=["bus", "carrier"])
    .droplevel("component")
    .loc[REGIONS]
    .drop("HVDC", level="carrier")
)

link_flow = n.links_t.p0.mean()[n.links.carrier == "HVDC"]

link_width = n.links.p_nom_opt[n.links.carrier == "HVDC"]

link_loading = n.links_t.p0.abs().mean()[n.links.carrier == "HVDC"] / link_width

n.plot(
    bus_size=bus_size / 4e7,
    margin=0.75,
    bus_split_circle=True,
    link_flow=link_flow / 50,
    link_width=link_width / 1e3,
    link_color=link_loading,
)

bus_size = ( n.statistics.energy_balance(bus_carrier="AC", groupby=["bus", "carrier"]) .droplevel("component") .loc[REGIONS] .drop("HVDC", level="carrier") ) link_flow = n.links_t.p0.mean()[n.links.carrier == "HVDC"] link_width = n.links.p_nom_opt[n.links.carrier == "HVDC"] link_loading = n.links_t.p0.abs().mean()[n.links.carrier == "HVDC"] / link_width n.plot( bus_size=bus_size / 4e7, margin=0.75, bus_split_circle=True, link_flow=link_flow / 50, link_width=link_width / 1e3, link_color=link_loading, )

Out[20]:

{'nodes': {'Bus': <matplotlib.collections.PatchCollection at 0x7a0b13b7b250>},
 'branches': {'Link': <matplotlib.collections.LineCollection at 0x7a0b13b7b390>},
 'flows': {'Link': <matplotlib.collections.PatchCollection at 0x7a0b13b7b610>}}

Of course, we can also just retrieve the optimised interconnector capacities as pandas.Series.

In [21]:

n.links.query("carrier == 'HVDC'").p_nom_opt.round(2)

n.links.query("carrier == 'HVDC'").p_nom_opt.round(2)

Out[21]:

name
VIC-NSW    2512.60
VIC-SA     2349.64
NSW-SA      434.72
Name: p_nom_opt, dtype: float64

A further interesting insight reveals itself when we plto the interconnector flows over the year (here as weekly rolling averages). The seasonality of the link between Victoria and New South Wales is striking, with VIC exporting to NSW in winter and importing in summer (Southern Hemisphere!).

In [22]:

n.links_t.p0.loc[:, n.links.carrier == "HVDC"].rolling("7d").mean().plot(
    ylim=[-2500, 2500], grid=True
)

n.links_t.p0.loc[:, n.links.carrier == "HVDC"].rolling("7d").mean().plot( ylim=[-2500, 2500], grid=True )

Out[22]:

<Axes: xlabel='snapshot'>

Also the state of charge profile of the hydrogen storage shows a clear seasonal pattern, with the storage being full in summer and empty in winter (Southern Hemisphere!).

In [23]:

n.stores_t.e.filter(like="hydrogen").sum(axis=1).div(1e6).plot(
    ylabel="Hydrogen storage [TWh]", grid=True
)

n.stores_t.e.filter(like="hydrogen").sum(axis=1).div(1e6).plot( ylabel="Hydrogen storage [TWh]", grid=True )

Out[23]:

<Axes: xlabel='snapshot', ylabel='Hydrogen storage [TWh]'>

Finally, there are some more functions to interactively plot hourly energy balance time series for different carriers, such as electricity:

In [24]:

n.statistics.energy_balance.iplot.area(bus_carrier="AC")

n.statistics.energy_balance.iplot.area(bus_carrier="AC")