In [1]:

import warnings

warnings.filterwarnings(
    "ignore"
)  # For better visibility, warnings are turned off in this notebook

import warnings warnings.filterwarnings( "ignore" ) # For better visibility, warnings are turned off in this notebook

ETHOS.FINE Tutorial: 2-nodal Electricity Supply System¶

In this application of the ETHOS.FINE framework, an energy supply system, consisting of two-regions, is modeled and optimized.

The workflow is structures as follows:

• Required packages are imported
• An energy system model instance is created
• Commodity sources are added to the energy supply system model
• Commodity conversion components are added to the energy supply system model
• Commodity storages are added to the energy supply system model
• Commodity transmission components are added to the energy supply system model
• Commodity sinks are added to the energy supply system model
• The energy supply system model is optimized
• Selected optimization results are presented

Import required packages¶

The ETHOS.FINE framework is imported which provides the required classes and functions for modeling the energy system.

In [2]:

import fine as fn  # Provides objects and functions to model an energy system
import pandas as pd  # Used to manage data in tables
import shapely as shp  # Used to generate geometric objects
import numpy as np  # Used to generate random input data

np.random.seed(
    42
)  # Sets a "seed" to produce the same random input data in each model run

import fine as fn # Provides objects and functions to model an energy system import pandas as pd # Used to manage data in tables import shapely as shp # Used to generate geometric objects import numpy as np # Used to generate random input data np.random.seed( 42 ) # Sets a "seed" to produce the same random input data in each model run

In [3]:

nbval-skip

import geopandas as gpd  # Used to display geo-referenced plots

import geopandas as gpd # Used to display geo-referenced plots

Model an energy system¶

Create an energy system model instance¶

The structure of the energy supply system model is given by the considered locations, commodities, the number of time steps as well as the hours per time step.

The commodities are specified by a unit (i.e. 'GW_electric', 'GW_naturalGas_lowerHeatingValue', 'Mio. t CO2/h') which can be given as an energy or mass unit per hour. Furthermore, the cost unit and length unit are specified.

In [4]:

# Input parameters
locations = {"regionN", "regionS"}
commodityUnitDict = {
    "electricity": r"GW$_{el}$",
    "naturalGas": r"GW$_{CH_{4},LHV}$",
    "CO2": r"Mio. t$_{CO_2}$/h",
}
commodities = {"electricity", "naturalGas", "CO2"}
numberOfTimeSteps, hoursPerTimeStep = 8760, 1
costUnit, lengthUnit = "1e6 Euro", "km"
CO2_reductionTarget = 0.8
balanceLimit = pd.DataFrame(
    index=["CO2 limit"],
    columns=["Total", "lowerBound"],
    data=[[-366 * (1 - CO2_reductionTarget), True]],
)

# Code
esM = fn.EnergySystemModel(
    locations=locations,
    commodities=commodities,
    numberOfTimeSteps=numberOfTimeSteps,
    commodityUnitsDict=commodityUnitDict,
    hoursPerTimeStep=hoursPerTimeStep,
    costUnit=costUnit,
    lengthUnit=lengthUnit,
    verboseLogLevel=0,
    balanceLimit=balanceLimit,
)

# Input parameters locations = {"regionN", "regionS"} commodityUnitDict = { "electricity": r"GW$_{el}$", "naturalGas": r"GW$_{CH_{4},LHV}$", "CO2": r"Mio. t$_{CO_2}$/h", } commodities = {"electricity", "naturalGas", "CO2"} numberOfTimeSteps, hoursPerTimeStep = 8760, 1 costUnit, lengthUnit = "1e6 Euro", "km" CO2_reductionTarget = 0.8 balanceLimit = pd.DataFrame( index=["CO2 limit"], columns=["Total", "lowerBound"], data=[[-366 * (1 - CO2_reductionTarget), True]], ) # Code esM = fn.EnergySystemModel( locations=locations, commodities=commodities, numberOfTimeSteps=numberOfTimeSteps, commodityUnitsDict=commodityUnitDict, hoursPerTimeStep=hoursPerTimeStep, costUnit=costUnit, lengthUnit=lengthUnit, verboseLogLevel=0, balanceLimit=balanceLimit, )

Add source components¶

Source components generate commodities across the energy system's virtual boundaries.

In [5]:

# Input parameters
name, commodity = "Wind turbines", "electricity"
hasCapacityVariable = True
operationRateMax = pd.DataFrame(
    [[np.random.beta(a=2, b=7.5), np.random.beta(a=2, b=9)] for t in range(8760)],
    index=range(8760),
    columns=["regionN", "regionS"],
).round(6)
capacityMax = pd.Series([400, 200], index=["regionN", "regionS"])
investPerCapacity, opexPerCapacity = 1200, 1200 * 0.02
interestRate, economicLifetime = 0.08, 20

# If data should be read from an excel file:
writer = pd.ExcelWriter("windTurbineProfile.xlsx")  # writes data to an excel file
operationRateMax.to_excel(writer)  # (not required if excel file
writer.close()  #  already exists)
operationRateMax = pd.read_excel(
    "windTurbineProfile.xlsx", index_col=0
)  # reads an excel file located in
# the current working directory

# Code
esM.add(
    fn.Source(
        esM=esM,
        name=name,
        commodity=commodity,
        hasCapacityVariable=hasCapacityVariable,
        operationRateMax=operationRateMax,
        capacityMax=capacityMax,
        investPerCapacity=investPerCapacity,
        opexPerCapacity=opexPerCapacity,
        interestRate=interestRate,
        economicLifetime=economicLifetime,
    )
)

# Input parameters name, commodity = "Wind turbines", "electricity" hasCapacityVariable = True operationRateMax = pd.DataFrame( [[np.random.beta(a=2, b=7.5), np.random.beta(a=2, b=9)] for t in range(8760)], index=range(8760), columns=["regionN", "regionS"], ).round(6) capacityMax = pd.Series([400, 200], index=["regionN", "regionS"]) investPerCapacity, opexPerCapacity = 1200, 1200 * 0.02 interestRate, economicLifetime = 0.08, 20 # If data should be read from an excel file: writer = pd.ExcelWriter("windTurbineProfile.xlsx") # writes data to an excel file operationRateMax.to_excel(writer) # (not required if excel file writer.close() # already exists) operationRateMax = pd.read_excel( "windTurbineProfile.xlsx", index_col=0 ) # reads an excel file located in # the current working directory # Code esM.add( fn.Source( esM=esM, name=name, commodity=commodity, hasCapacityVariable=hasCapacityVariable, operationRateMax=operationRateMax, capacityMax=capacityMax, investPerCapacity=investPerCapacity, opexPerCapacity=opexPerCapacity, interestRate=interestRate, economicLifetime=economicLifetime, ) )

In [6]:

# Input parameters
name, commodity = "PV", "electricity"
hasCapacityVariable = True
dailyProfileSimple = [
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0.05,
    0.15,
    0.2,
    0.4,
    0.8,
    0.7,
    0.4,
    0.2,
    0.15,
    0.05,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
]
operationRateMax = pd.DataFrame(
    [[u, u] for day in range(365) for u in dailyProfileSimple],
    index=range(8760),
    columns=["regionN", "regionS"],
)
capacityMax = pd.Series([100, 100], index=["regionN", "regionS"])
investPerCapacity, opexPerCapacity = 800, 800 * 0.02
interestRate, economicLifetime = 0.08, 25

# If data should be read from an excel file:
writer = pd.ExcelWriter("PV_Profile.xlsx")  # writes data to an excel file
operationRateMax.to_excel(writer)  # (not required if excel file
writer.close()  #  already exists)
operationRateMax = pd.read_excel(
    "PV_Profile.xlsx", index_col=0
)  # reads an excel file located in
# the current working directory

# Code
esM.add(
    fn.Source(
        esM=esM,
        name=name,
        commodity=commodity,
        hasCapacityVariable=hasCapacityVariable,
        operationRateMax=operationRateMax,
        capacityMax=capacityMax,
        investPerCapacity=investPerCapacity,
        opexPerCapacity=opexPerCapacity,
        interestRate=interestRate,
        economicLifetime=economicLifetime,
    )
)

# Input parameters name, commodity = "PV", "electricity" hasCapacityVariable = True dailyProfileSimple = [ 0, 0, 0, 0, 0, 0, 0, 0.05, 0.15, 0.2, 0.4, 0.8, 0.7, 0.4, 0.2, 0.15, 0.05, 0, 0, 0, 0, 0, 0, 0, ] operationRateMax = pd.DataFrame( [[u, u] for day in range(365) for u in dailyProfileSimple], index=range(8760), columns=["regionN", "regionS"], ) capacityMax = pd.Series([100, 100], index=["regionN", "regionS"]) investPerCapacity, opexPerCapacity = 800, 800 * 0.02 interestRate, economicLifetime = 0.08, 25 # If data should be read from an excel file: writer = pd.ExcelWriter("PV_Profile.xlsx") # writes data to an excel file operationRateMax.to_excel(writer) # (not required if excel file writer.close() # already exists) operationRateMax = pd.read_excel( "PV_Profile.xlsx", index_col=0 ) # reads an excel file located in # the current working directory # Code esM.add( fn.Source( esM=esM, name=name, commodity=commodity, hasCapacityVariable=hasCapacityVariable, operationRateMax=operationRateMax, capacityMax=capacityMax, investPerCapacity=investPerCapacity, opexPerCapacity=opexPerCapacity, interestRate=interestRate, economicLifetime=economicLifetime, ) )

In [7]:

# Input parameters
name, commodity = "Natural gas import", "naturalGas"
hasCapacityVariable = False
commodityCost = 0.03

# Code
esM.add(
    fn.Source(
        esM=esM,
        name=name,
        commodity=commodity,
        hasCapacityVariable=hasCapacityVariable,
        commodityCost=commodityCost,
    )
)

# Input parameters name, commodity = "Natural gas import", "naturalGas" hasCapacityVariable = False commodityCost = 0.03 # Code esM.add( fn.Source( esM=esM, name=name, commodity=commodity, hasCapacityVariable=hasCapacityVariable, commodityCost=commodityCost, ) )

Add conversion components¶

Conversion components convert m commodities into n other commodities.

In [8]:

# Input parameters
name, physicalUnit = "Gas power plants", r"GW$_{el}$"
commodityConversionFactors = {
    "electricity": 1,
    "naturalGas": -1 / 0.63,
    "CO2": 201 * 1e-6 / 0.63,
}
hasCapacityVariable = True
investPerCapacity, opexPerCapacity = 650, 650 * 0.03
interestRate, economicLifetime = 0.08, 30

# Code
esM.add(
    fn.Conversion(
        esM=esM,
        name=name,
        physicalUnit=physicalUnit,
        commodityConversionFactors=commodityConversionFactors,
        hasCapacityVariable=hasCapacityVariable,
        investPerCapacity=investPerCapacity,
        opexPerCapacity=opexPerCapacity,
        interestRate=interestRate,
        economicLifetime=economicLifetime,
    )
)

# Input parameters name, physicalUnit = "Gas power plants", r"GW$_{el}$" commodityConversionFactors = { "electricity": 1, "naturalGas": -1 / 0.63, "CO2": 201 * 1e-6 / 0.63, } hasCapacityVariable = True investPerCapacity, opexPerCapacity = 650, 650 * 0.03 interestRate, economicLifetime = 0.08, 30 # Code esM.add( fn.Conversion( esM=esM, name=name, physicalUnit=physicalUnit, commodityConversionFactors=commodityConversionFactors, hasCapacityVariable=hasCapacityVariable, investPerCapacity=investPerCapacity, opexPerCapacity=opexPerCapacity, interestRate=interestRate, economicLifetime=economicLifetime, ) )

Add storage components¶

Storage components can store commodities across time steps.

The self discharge of a storage technology is described in FINE in percent per hour. If the literature value is given in percent per month, e.g. 3%/month, the self discharge per hours is obtained using the equation (1-$\text{selfDischarge}_\text{hour})^{30*24\text{h}} = 1-\text{selfDischarge}_\text{month}$.

In [9]:

# Input parameters
name, commodity = "Batteries", "electricity"
hasCapacityVariable = True
chargeEfficiency, dischargeEfficiency, selfDischarge = (
    0.95,
    0.95,
    1 - (1 - 0.03) ** (1 / (30 * 24)),
)
chargeRate, dischargeRate = 1, 1
investPerCapacity, opexPerCapacity = 150, 150 * 0.01
interestRate, economicLifetime, cyclicLifetime = 0.08, 22, 12000

# Code
esM.add(
    fn.Storage(
        esM=esM,
        name=name,
        commodity=commodity,
        hasCapacityVariable=hasCapacityVariable,
        chargeEfficiency=chargeEfficiency,
        cyclicLifetime=cyclicLifetime,
        dischargeEfficiency=dischargeEfficiency,
        selfDischarge=selfDischarge,
        chargeRate=chargeRate,
        dischargeRate=dischargeRate,
        investPerCapacity=investPerCapacity,
        opexPerCapacity=opexPerCapacity,
        interestRate=interestRate,
        economicLifetime=economicLifetime,
    )
)

# Input parameters name, commodity = "Batteries", "electricity" hasCapacityVariable = True chargeEfficiency, dischargeEfficiency, selfDischarge = ( 0.95, 0.95, 1 - (1 - 0.03) ** (1 / (30 * 24)), ) chargeRate, dischargeRate = 1, 1 investPerCapacity, opexPerCapacity = 150, 150 * 0.01 interestRate, economicLifetime, cyclicLifetime = 0.08, 22, 12000 # Code esM.add( fn.Storage( esM=esM, name=name, commodity=commodity, hasCapacityVariable=hasCapacityVariable, chargeEfficiency=chargeEfficiency, cyclicLifetime=cyclicLifetime, dischargeEfficiency=dischargeEfficiency, selfDischarge=selfDischarge, chargeRate=chargeRate, dischargeRate=dischargeRate, investPerCapacity=investPerCapacity, opexPerCapacity=opexPerCapacity, interestRate=interestRate, economicLifetime=economicLifetime, ) )

Add transmission components¶

Transmission components transmit commodities between regions.

In [10]:

# Input parameters
name, commodity = "AC cables", "electricity"
hasCapacityVariable = True
capacityFix = pd.DataFrame(
    [[0, 30], [30, 0]], columns=["regionN", "regionS"], index=["regionN", "regionS"]
)
distances = pd.DataFrame(
    [[0, 400], [400, 0]], columns=["regionN", "regionS"], index=["regionN", "regionS"]
)
losses = 0.0001

# Code
esM.add(
    fn.Transmission(
        esM=esM,
        name=name,
        commodity=commodity,
        hasCapacityVariable=hasCapacityVariable,
        capacityFix=capacityFix,
        distances=distances,
        losses=losses,
    )
)

# Input parameters name, commodity = "AC cables", "electricity" hasCapacityVariable = True capacityFix = pd.DataFrame( [[0, 30], [30, 0]], columns=["regionN", "regionS"], index=["regionN", "regionS"] ) distances = pd.DataFrame( [[0, 400], [400, 0]], columns=["regionN", "regionS"], index=["regionN", "regionS"] ) losses = 0.0001 # Code esM.add( fn.Transmission( esM=esM, name=name, commodity=commodity, hasCapacityVariable=hasCapacityVariable, capacityFix=capacityFix, distances=distances, losses=losses, ) )

In [11]:

distances

distances

Out[11]:

	regionN	regionS
regionN	0	400
regionS	400	0

Add sink components¶

Sinks remove commodities across the energy system´s virtual boundaries.

In [12]:

# Input parameters
name, commodity = (
    "Electricity demand",
    "electricity",
)
hasCapacityVariable = False
dailyProfileSimple = [
    0.6,
    0.6,
    0.6,
    0.6,
    0.6,
    0.7,
    0.9,
    1,
    1,
    1,
    1,
    1,
    1,
    1,
    1,
    1,
    1,
    1,
    1,
    1,
    1,
    1,
    0.9,
    0.8,
]
operationRateFix = pd.DataFrame(
    [
        [(u + 0.1 * np.random.rand()) * 25, (u + 0.1 * np.random.rand()) * 40]
        for day in range(365)
        for u in dailyProfileSimple
    ],
    index=range(8760),
    columns=["regionN", "regionS"],
).round(2)

# If data should be read from an excel file:
writer = pd.ExcelWriter("demandProfile.xlsx")  # writes data to an excel file
operationRateFix.to_excel(writer)  # (not required if excel file
writer.close()  #  already exists)
operationRateFix = pd.read_excel(
    "demandProfile.xlsx", index_col=0
)  # reads an excel file located in
# the current working directory

# Code
esM.add(
    fn.Sink(
        esM=esM,
        name=name,
        commodity=commodity,
        hasCapacityVariable=hasCapacityVariable,
        operationRateFix=operationRateFix,
    )
)

# Input parameters name, commodity = ( "Electricity demand", "electricity", ) hasCapacityVariable = False dailyProfileSimple = [ 0.6, 0.6, 0.6, 0.6, 0.6, 0.7, 0.9, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0.9, 0.8, ] operationRateFix = pd.DataFrame( [ [(u + 0.1 * np.random.rand()) * 25, (u + 0.1 * np.random.rand()) * 40] for day in range(365) for u in dailyProfileSimple ], index=range(8760), columns=["regionN", "regionS"], ).round(2) # If data should be read from an excel file: writer = pd.ExcelWriter("demandProfile.xlsx") # writes data to an excel file operationRateFix.to_excel(writer) # (not required if excel file writer.close() # already exists) operationRateFix = pd.read_excel( "demandProfile.xlsx", index_col=0 ) # reads an excel file located in # the current working directory # Code esM.add( fn.Sink( esM=esM, name=name, commodity=commodity, hasCapacityVariable=hasCapacityVariable, operationRateFix=operationRateFix, ) )

In [13]:

# Input parameters
name, commodity = (
    "CO2 to enviroment",
    "CO2",
)
hasCapacityVariable = False

# Code
esM.add(
    fn.Sink(
        esM=esM,
        name=name,
        commodity=commodity,
        hasCapacityVariable=hasCapacityVariable,
        balanceLimitID="CO2 limit",
    )
)

# Input parameters name, commodity = ( "CO2 to enviroment", "CO2", ) hasCapacityVariable = False # Code esM.add( fn.Sink( esM=esM, name=name, commodity=commodity, hasCapacityVariable=hasCapacityVariable, balanceLimitID="CO2 limit", ) )

Optimize energy system model¶

All components are now added to the model and the model can be optimized. If the computational complexity of the optimization should be reduced, the time series data of the specified components can be clustered before the optimization and the parameter timeSeriesAggregation is set to True in the optimize call.

In [14]:

# Input parameters
numberOfTypicalPeriods = 30

# Code
esM.aggregateTemporally(numberOfTypicalPeriods=numberOfTypicalPeriods)

# Input parameters numberOfTypicalPeriods = 30 # Code esM.aggregateTemporally(numberOfTypicalPeriods=numberOfTypicalPeriods)

Clustering time series data with 30 typical periods and 24 time steps per period 
further clustered to 12 segments per period...
		(5.3689 sec)

In [15]:

# Input parameters
timeSeriesAggregation = True
solver = fn.utils.ImplementedSolvers.STANDARD_SOLVER.value

# Code
esM.optimize(timeSeriesAggregation=timeSeriesAggregation, solver=solver)

# Input parameters timeSeriesAggregation = True solver = fn.utils.ImplementedSolvers.STANDARD_SOLVER.value # Code esM.optimize(timeSeriesAggregation=timeSeriesAggregation, solver=solver)

Time series aggregation specifications:
Number of typical periods:30, number of time steps per period:24, number of segments per period:12

Declaring sets, variables and constraints for SourceSinkModel
	declaring sets... 
	declaring variables... 
	declaring constraints... 
		(0.5833 sec)

Declaring sets, variables and constraints for ConversionModel
	declaring sets... 
	declaring variables... 
	declaring constraints... 
		(0.0964 sec)

Declaring sets, variables and constraints for StorageModel
	declaring sets... 
	declaring variables... 
	declaring constraints... 
		(0.3419 sec)

Declaring sets, variables and constraints for TransmissionModel
	declaring sets... 
	declaring variables... 
	declaring constraints... 
		(0.0530 sec)

Declaring shared potential constraint...
		(0.0001 sec)

Declaring linked component quantity constraint...
		(0.0000 sec)

Declaring commodity balances...
		(0.0527 sec)

		(0.0020 sec)

Declaring objective function...
		(0.5272 sec)

GLPSOL--GLPK LP/MIP Solver 5.0
Parameter(s) specified in the command line:
 --write /tmp/tmpy2qvi_0b.glpk.raw --wglp /tmp/tmpnmn5dm3b.glpk.glp --cpxlp
 /tmp/tmpdo_yrc2_.pyomo.lp
Reading problem data from '/tmp/tmpdo_yrc2_.pyomo.lp'...
11599 rows, 8152 columns, 32238 non-zeros
75924 lines were read
Writing problem data to '/tmp/tmpnmn5dm3b.glpk.glp'...
61571 lines were written
GLPK Simplex Optimizer 5.0
11599 rows, 8152 columns, 32238 non-zeros
Preprocessing...
10391 rows, 6286 columns, 28324 non-zeros
Scaling...
 A: min|aij| =  3.190e-04  max|aij| =  5.455e+02  ratio =  1.710e+06
GM: min|aij| =  2.926e-01  max|aij| =  3.417e+00  ratio =  1.168e+01
EQ: min|aij| =  8.563e-02  max|aij| =  1.000e+00  ratio =  1.168e+01
Constructing initial basis...
Size of triangular part is 9659
      0: obj =   0.000000000e+00 inf =   1.390e+04 (720)
    743: obj =   9.946047517e+04 inf =   0.000e+00 (0) 7
Perturbing LP to avoid instability [1688]...
Removing LP perturbation [6161]...
*  6161: obj =   3.976374573e+04 inf =   8.052e-12 (0) 49
OPTIMAL LP SOLUTION FOUND
Time used:   3.1 secs
Memory used: 14.7 Mb (15457837 bytes)
Writing basic solution to '/tmp/tmpy2qvi_0b.glpk.raw'...
19760 lines were written

Status: ok
Termination condition: optimal
Statistics: 
  Branch and bound: 
    Number of bounded subproblems: 0
    Number of created subproblems: 0
Error rc: 0
Time: 3.1877660751342773


Name: unknown
Lower bound: 39763.7457348323
Upper bound: 39763.7457348323
Number of objectives: 1
Number of constraints: 11599
Number of variables: 8152
Number of nonzeros: 32238
Sense: minimize

Solve time: 3.4588963985443115 sec.

Processing optimization output...
for SourceSinkModel ...    (0.8489sec)
for ConversionModel ...    (0.2056sec)
for StorageModel ...       (0.6905sec)
for TransmissionModel ...  (0.3261sec)
		(2.0907 sec)

Output of selected results¶

For the assessment of the optimization result, several result output functions are available. They can be categorized into output in form of tables, geo-referenced output visualization and the full time series visualization.

Create a regional shape file and visualize it¶

Information on the geometrical shape of the investigated regions can either be downloaded from a website (e.g. from https://gadm.org/) or manually created. In this notebook, the geometries are manually created.

In [16]:

nbval-skip

# Create two circles, representing the two regions, and store their geometries in a shape file
shpRegionS = shp.geometry.Point(0.5, 0.5).buffer(0.5)
shpRegionN = shp.geometry.Point(0.5, 1.5).buffer(0.5)
regionsGdf = gpd.GeoDataFrame(
    {"geometry": [shpRegionN, shpRegionS], "regionName": ["regionN", "regionS"]},
    index=["regionN", "regionS"],
    crs="epsg:3035",
)
regionsGdf.to_file("regions.shp")

# Create a line, representing the connection between the two regions, and store its geometry in a
# shape file
lines = shp.geometry.LineString([(0.5, 0.5), (0.5, 1.5)])
linesGdf = gpd.GeoDataFrame(
    {
        "geometry": [lines, lines],
        "loc0": ["regionN", "regionS"],
        "loc1": ["regionS", "regionN"],
    },
    index=["regionN_regionS", "regionS_regionN"],
    crs="epsg:3035",
)
linesGdf.to_file("lines.shp")

# Visualize the geometric representation of the two regions
fig, ax = fn.plotLocations("regions.shp", indexColumn="regionName", plotLocNames=True)

# Create two circles, representing the two regions, and store their geometries in a shape file shpRegionS = shp.geometry.Point(0.5, 0.5).buffer(0.5) shpRegionN = shp.geometry.Point(0.5, 1.5).buffer(0.5) regionsGdf = gpd.GeoDataFrame( {"geometry": [shpRegionN, shpRegionS], "regionName": ["regionN", "regionS"]}, index=["regionN", "regionS"], crs="epsg:3035", ) regionsGdf.to_file("regions.shp") # Create a line, representing the connection between the two regions, and store its geometry in a # shape file lines = shp.geometry.LineString([(0.5, 0.5), (0.5, 1.5)]) linesGdf = gpd.GeoDataFrame( { "geometry": [lines, lines], "loc0": ["regionN", "regionS"], "loc1": ["regionS", "regionN"], }, index=["regionN_regionS", "regionS_regionN"], crs="epsg:3035", ) linesGdf.to_file("lines.shp") # Visualize the geometric representation of the two regions fig, ax = fn.plotLocations("regions.shp", indexColumn="regionName", plotLocNames=True)

Display optimization summaries¶

For each modeling class, an optimization summary can be stored and displayed.

In [17]:

srcSnkSummary = esM.getOptimizationSummary("SourceSinkModel", outputLevel=1)
display(esM.getOptimizationSummary("SourceSinkModel", outputLevel=2))

srcSnkSummary = esM.getOptimizationSummary("SourceSinkModel", outputLevel=1) display(esM.getOptimizationSummary("SourceSinkModel", outputLevel=2))

			regionN	regionS
Component	Property	Unit
CO2 to enviroment	operation	[Mio. t$_{CO_2}$/h*h/a]	0.517539	72.682461
		[Mio. t$_{CO_2}$/h*h]	0.517539	72.682461
Electricity demand	operation	[GW$_{el}$*h/a]	205207.79	328596.33
		[GW$_{el}$*h]	205207.79	328596.33
Natural gas import	NPVcontribution	[1e6 Euro]	77.244598	10848.128536
	TAC	[1e6 Euro/a]	77.244598	10848.128536
	commodCosts	[1e6 Euro/a]	77.244598	10848.128536
	operation	[GW$_{CH_{4},LHV}$*h/a]	2574.81994	361604.284537
		[GW$_{CH_{4},LHV}$*h]	2574.81994	361604.284537
PV	NPVcontribution	[1e6 Euro]	0.0	1959.001998
	TAC	[1e6 Euro/a]	0.0	1959.001998
	capacity	[GW$_{el}$]	0.0	21.540982
	capexCap	[1e6 Euro/a]	0.0	1614.346291
	commissioning	[GW$_{el}$]	-0.0	21.540982
	invest	[1e6 Euro]	0.0	17232.785352
	operation	[GW$_{el}$*h/a]	0.0	24373.620782
		[GW$_{el}$*h]	0.0	24373.620782
	opexCap	[1e6 Euro/a]	0.0	344.655707
Wind turbines	NPVcontribution	[1e6 Euro]	22638.368728	0.0
	TAC	[1e6 Euro/a]	22638.368728	0.0
	capacity	[GW$_{el}$]	154.821217	0.0
	capexCap	[1e6 Euro/a]	18922.659518	0.0
	commissioning	[GW$_{el}$]	154.821217	-0.0
	invest	[1e6 Euro]	185785.460493	0.0
	operation	[GW$_{el}$*h/a]	285458.470678	0.0
		[GW$_{el}$*h]	285458.470678	0.0
	opexCap	[1e6 Euro/a]	3715.70921	0.0

In [18]:

convSummary = esM.getOptimizationSummary("ConversionModel", outputLevel=1)
display(esM.getOptimizationSummary("ConversionModel", outputLevel=2))

convSummary = esM.getOptimizationSummary("ConversionModel", outputLevel=1) display(esM.getOptimizationSummary("ConversionModel", outputLevel=2))

			regionN	regionS
Component	Property	Unit
Gas power plants	NPVcontribution	[1e6 Euro]	51.482796	3138.932574
	TAC	[1e6 Euro/a]	51.482796	3138.932574
	capacity	[GW$_{el}$]	0.666549	40.639833
	capexCap	[1e6 Euro/a]	38.485091	2346.455833
	commissioning	[GW$_{el}$]	0.666549	40.639833
	invest	[1e6 Euro]	433.256822	26415.891387
	operation	[GW$_{el}$*h/a]	1622.136563	227810.699258
		[GW$_{el}$*h]	1622.136563	227810.699258
	opexCap	[1e6 Euro/a]	12.997705	792.476742

In [19]:

storSummary = esM.getOptimizationSummary("StorageModel", outputLevel=1)
display(esM.getOptimizationSummary("StorageModel", outputLevel=2))

storSummary = esM.getOptimizationSummary("StorageModel", outputLevel=1) display(esM.getOptimizationSummary("StorageModel", outputLevel=2))

			regionN	regionS
Component	Property	Unit
Batteries	NPVcontribution	[1e6 Euro]	979.809612	70.776892
	TAC	[1e6 Euro/a]	979.809612	70.776892
	capacity	[GW$_{el}$*h]	60.464121	4.367647
	capexCap	[1e6 Euro/a]	889.11343	64.225421
	commissioning	[GW$_{el}$*h]	60.464121	4.367647
	invest	[1e6 Euro]	9069.618184	655.147059
	operationCharge	[GW$_{el}$*h/a]	18027.11586	425.732468
		[GW$_{el}$*h]	18027.11586	425.732468
	operationDischarge	[GW$_{el}$*h/a]	16259.862567	384.040369
		[GW$_{el}$*h]	16259.862567	384.040369
	opexCap	[1e6 Euro/a]	90.696182	6.551471

In [20]:

transSummary = esM.getOptimizationSummary("TransmissionModel", outputLevel=1)
display(esM.getOptimizationSummary("TransmissionModel", outputLevel=2))

transSummary = esM.getOptimizationSummary("TransmissionModel", outputLevel=1) display(esM.getOptimizationSummary("TransmissionModel", outputLevel=2))

				regionN	regionS
Component	Property	Unit	LocationIn
AC cables	capacity	[GW$_{el}$]	regionN	NaN	30.0
			regionS	30.0	NaN
	commissioning	[GW$_{el}$]	regionN	NaN	30.0
			regionS	30.0	NaN
	operation	[GW$_{el}$*h/a]	regionN	NaN	85586.861879
			regionS	5709.685346	NaN
		[GW$_{el}$*h]	regionN	NaN	85586.861879
			regionS	5709.685346	NaN

Display regional and operational plots¶

Georeferenced plots as well as plots representing time series can be displayed for each component.

Wind turbines¶

In [21]:

nbval-skip

# If wind turbines are built, their capacities are displayed in a geo-referenced plot
if srcSnkSummary.loc[("Wind turbines", "capacity", "[GW$_{el}$]")].sum() > 0:
    fig, ax = fn.plotLocationalColorMap(
        esM,
        "Wind turbines",
        "regions.shp",
        "regionName",
        perArea=False,
        zlabel="Capacity\n[GW]\n",
        figsize=(4, 4),
    )
else:
    print("No wind turbines built.")

# If wind turbines are built in regionN, their operation is displayed as heatmap
if srcSnkSummary.loc[("Wind turbines", "capacity", "[GW$_{el}$]"), "regionN"] > 0:
    fig, ax = fn.plotOperationColorMap(
        esM,
        "Wind turbines",
        "regionN",
        figsize=(5, 3),
        xlabel="Day of the year",
        ylabel="Hour of the day",
        zlabel="Operation\nin regionN\n[GW]",
    )

# If wind turbines are built in regionS, their operation is displayed as heatmap
if srcSnkSummary.loc[("Wind turbines", "capacity", "[GW$_{el}$]"), "regionS"] > 0:
    fig, ax = fn.plotOperationColorMap(
        esM,
        "Wind turbines",
        "regionS",
        figsize=(5, 3),
        xlabel="Day of the year",
        ylabel="Hour of the Day",
        zlabel="Operation\nin regionS\n[GW]",
        orientation="vertical",
    )

# If wind turbines are built, their capacities are displayed in a geo-referenced plot if srcSnkSummary.loc[("Wind turbines", "capacity", "[GW$_{el}$]")].sum() > 0: fig, ax = fn.plotLocationalColorMap( esM, "Wind turbines", "regions.shp", "regionName", perArea=False, zlabel="Capacity\n[GW]\n", figsize=(4, 4), ) else: print("No wind turbines built.") # If wind turbines are built in regionN, their operation is displayed as heatmap if srcSnkSummary.loc[("Wind turbines", "capacity", "[GW$_{el}$]"), "regionN"] > 0: fig, ax = fn.plotOperationColorMap( esM, "Wind turbines", "regionN", figsize=(5, 3), xlabel="Day of the year", ylabel="Hour of the day", zlabel="Operation\nin regionN\n[GW]", ) # If wind turbines are built in regionS, their operation is displayed as heatmap if srcSnkSummary.loc[("Wind turbines", "capacity", "[GW$_{el}$]"), "regionS"] > 0: fig, ax = fn.plotOperationColorMap( esM, "Wind turbines", "regionS", figsize=(5, 3), xlabel="Day of the year", ylabel="Hour of the Day", zlabel="Operation\nin regionS\n[GW]", orientation="vertical", )

PV systems¶

In [22]:

nbval-skip

# If PV systems are built, their capacities are displayed in a geo-referenced plot
if srcSnkSummary.loc[("PV", "capacity", "[GW$_{el}$]")].sum() > 0:
    fig, ax = fn.plotLocationalColorMap(
        esM,
        "PV",
        "regions.shp",
        "regionName",
        perArea=False,
        zlabel="Capacity\n[GW]\n",
        figsize=(4, 4),
    )
else:
    print("No PV systems built.")

# If PV systems are built in regionS, their operation is displayed as heatmap
if srcSnkSummary.loc[("PV", "capacity", "[GW$_{el}$]"), "regionN"] > 0:
    fig, ax = fn.plotOperationColorMap(
        esM,
        "PV",
        "regionN",
        figsize=(5, 3),
        xlabel="Day of the year",
        ylabel="Hour of the day",
        zlabel="Operation\nin regionN\n[GW]",
    )

# If PV systems are built in regionS, their operation is displayed as heatmap
if srcSnkSummary.loc[("PV", "capacity", "[GW$_{el}$]"), "regionS"] > 0:
    fig, ax = fn.plotOperationColorMap(
        esM,
        "PV",
        "regionS",
        figsize=(5, 3),
        xlabel="Day of the year",
        ylabel="Hour of the day",
        zlabel="Operation\nin regionS\n[GW]",
        orientation="vertical",
    )

# If PV systems are built, their capacities are displayed in a geo-referenced plot if srcSnkSummary.loc[("PV", "capacity", "[GW$_{el}$]")].sum() > 0: fig, ax = fn.plotLocationalColorMap( esM, "PV", "regions.shp", "regionName", perArea=False, zlabel="Capacity\n[GW]\n", figsize=(4, 4), ) else: print("No PV systems built.") # If PV systems are built in regionS, their operation is displayed as heatmap if srcSnkSummary.loc[("PV", "capacity", "[GW$_{el}$]"), "regionN"] > 0: fig, ax = fn.plotOperationColorMap( esM, "PV", "regionN", figsize=(5, 3), xlabel="Day of the year", ylabel="Hour of the day", zlabel="Operation\nin regionN\n[GW]", ) # If PV systems are built in regionS, their operation is displayed as heatmap if srcSnkSummary.loc[("PV", "capacity", "[GW$_{el}$]"), "regionS"] > 0: fig, ax = fn.plotOperationColorMap( esM, "PV", "regionS", figsize=(5, 3), xlabel="Day of the year", ylabel="Hour of the day", zlabel="Operation\nin regionS\n[GW]", orientation="vertical", )

Gas power plants¶

In [23]:

nbval-skip

# If CCGT plants are built, their capacities are displayed in a geo-referenced plot
if convSummary.loc[("Gas power plants", "capacity", "[GW$_{el}$]")].sum() > 0:
    fig, ax = fn.plotLocationalColorMap(
        esM,
        "Gas power plants",
        "regions.shp",
        "regionName",
        perArea=False,
        zlabel="Capacity\n[GW]\n",
        figsize=(4, 4),
    )
else:
    print("No CCGT plants built.")

# If CCGT plants are built in regionS, their operation is displayed as heatmap
if convSummary.loc[("Gas power plants", "capacity", "[GW$_{el}$]"), "regionN"] > 0:
    fig, ax = fn.plotOperationColorMap(
        esM,
        "Gas power plants",
        "regionN",
        figsize=(5, 3),
        xlabel="Day of the year",
        ylabel="Hour of the day",
        zlabel="Operation\nin regionN\n[GW]",
        orientation="vertical",
    )

# If CCGT plants are built in regionS, their operation is displayed as heatmap
if convSummary.loc[("Gas power plants", "capacity", "[GW$_{el}$]"), "regionS"] > 0:
    fig, ax = fn.plotOperationColorMap(
        esM,
        "Gas power plants",
        "regionS",
        figsize=(5, 3),
        xlabel="Day of the year",
        ylabel="Hour of the day",
        zlabel="Operation\nin regionS\n[GW]",
        orientation="vertical",
    )

# If CCGT plants are built, their capacities are displayed in a geo-referenced plot if convSummary.loc[("Gas power plants", "capacity", "[GW$_{el}$]")].sum() > 0: fig, ax = fn.plotLocationalColorMap( esM, "Gas power plants", "regions.shp", "regionName", perArea=False, zlabel="Capacity\n[GW]\n", figsize=(4, 4), ) else: print("No CCGT plants built.") # If CCGT plants are built in regionS, their operation is displayed as heatmap if convSummary.loc[("Gas power plants", "capacity", "[GW$_{el}$]"), "regionN"] > 0: fig, ax = fn.plotOperationColorMap( esM, "Gas power plants", "regionN", figsize=(5, 3), xlabel="Day of the year", ylabel="Hour of the day", zlabel="Operation\nin regionN\n[GW]", orientation="vertical", ) # If CCGT plants are built in regionS, their operation is displayed as heatmap if convSummary.loc[("Gas power plants", "capacity", "[GW$_{el}$]"), "regionS"] > 0: fig, ax = fn.plotOperationColorMap( esM, "Gas power plants", "regionS", figsize=(5, 3), xlabel="Day of the year", ylabel="Hour of the day", zlabel="Operation\nin regionS\n[GW]", orientation="vertical", )

Batteries¶

In [24]:

nbval-skip

# If batteries are built, their capacities are displayed in a geo-referenced plot
if storSummary.loc[("Batteries", "capacity", "[GW$_{el}$*h]")].sum() > 0:
    fig, ax = fn.plotLocationalColorMap(
        esM,
        "Batteries",
        "regions.shp",
        "regionName",
        perArea=False,
        zlabel="Capacity\n[GWh]\n",
        figsize=(4, 4),
    )
else:
    print("No batteries built.")

# If batteries are built in regionS, their storage inventory is displayed as heatmap
if storSummary.loc[("Batteries", "capacity", "[GW$_{el}$*h]"), "regionN"] > 0:
    fig, ax = fn.plotOperationColorMap(
        esM,
        "Batteries",
        "regionN",
        figsize=(5, 3),
        xlabel="Day of the year",
        ylabel="Hour of the day",
        zlabel="State of charge\nin regionN\n[GW]",
        variableName="stateOfChargeOperationVariablesOptimum",
        orientation="vertical",
    )

# If batteries are built in regionS, their storage inventory is displayed as heatmap
if storSummary.loc[("Batteries", "capacity", "[GW$_{el}$*h]"), "regionS"] > 0:
    fig, ax = fn.plotOperationColorMap(
        esM,
        "Batteries",
        "regionS",
        figsize=(5, 3),
        xlabel="Day of the year",
        ylabel="Hour of the day",
        zlabel="State of charge\nin regionS\n[GW]",
        variableName="stateOfChargeOperationVariablesOptimum",
        orientation="vertical",
    )

# If batteries are built, their capacities are displayed in a geo-referenced plot if storSummary.loc[("Batteries", "capacity", "[GW$_{el}$*h]")].sum() > 0: fig, ax = fn.plotLocationalColorMap( esM, "Batteries", "regions.shp", "regionName", perArea=False, zlabel="Capacity\n[GWh]\n", figsize=(4, 4), ) else: print("No batteries built.") # If batteries are built in regionS, their storage inventory is displayed as heatmap if storSummary.loc[("Batteries", "capacity", "[GW$_{el}$*h]"), "regionN"] > 0: fig, ax = fn.plotOperationColorMap( esM, "Batteries", "regionN", figsize=(5, 3), xlabel="Day of the year", ylabel="Hour of the day", zlabel="State of charge\nin regionN\n[GW]", variableName="stateOfChargeOperationVariablesOptimum", orientation="vertical", ) # If batteries are built in regionS, their storage inventory is displayed as heatmap if storSummary.loc[("Batteries", "capacity", "[GW$_{el}$*h]"), "regionS"] > 0: fig, ax = fn.plotOperationColorMap( esM, "Batteries", "regionS", figsize=(5, 3), xlabel="Day of the year", ylabel="Hour of the day", zlabel="State of charge\nin regionS\n[GW]", variableName="stateOfChargeOperationVariablesOptimum", orientation="vertical", )

AC cables¶

In [25]:

nbval-skip

# The built AC cable capacities are displayed
fig, ax = fn.plotLocations("regions.shp", indexColumn="regionName")
fig, ax = fn.plotTransmission(
    esM,
    "AC cables",
    "lines.shp",
    loc0="loc0",
    loc1="loc1",
    fig=fig,
    ax=ax,
    cbHeight=0.4,
)

# The built AC cable capacities are displayed fig, ax = fn.plotLocations("regions.shp", indexColumn="regionName") fig, ax = fn.plotTransmission( esM, "AC cables", "lines.shp", loc0="loc0", loc1="loc1", fig=fig, ax=ax, cbHeight=0.4, )

Electricity demand¶

In [26]:

# The electricity demand time series in regionN is displayed
fig, ax = fn.plotOperationColorMap(
    esM,
    "Electricity demand",
    "regionN",
    figsize=(5, 3),
    xlabel="Day of the year",
    ylabel="Hour of the day",
    zlabel="Demand\nin regionN\n[GW]",
    orientation="vertical",
)

# The electricity demand time series in regionS is displayed
fig, ax = fn.plotOperationColorMap(
    esM,
    "Electricity demand",
    "regionS",
    figsize=(5, 3),
    xlabel="Day of the year",
    ylabel="Hour of the day",
    zlabel="Demand\nin regionS\n[GW]",
    orientation="vertical",
)

# The electricity demand time series in regionN is displayed fig, ax = fn.plotOperationColorMap( esM, "Electricity demand", "regionN", figsize=(5, 3), xlabel="Day of the year", ylabel="Hour of the day", zlabel="Demand\nin regionN\n[GW]", orientation="vertical", ) # The electricity demand time series in regionS is displayed fig, ax = fn.plotOperationColorMap( esM, "Electricity demand", "regionS", figsize=(5, 3), xlabel="Day of the year", ylabel="Hour of the day", zlabel="Demand\nin regionS\n[GW]", orientation="vertical", )

In [ ]: