Example notebook
[1]:
import matplotlib.pyplot as plt
from pommes import test_case_path
from pommes.io.build_input_dataset import (
build_input_parameters,
read_config_file,
)
from pommes.io.save_solution import save_solution
from pommes.model.build_model import build_model
from pommes.model.data_validation.dataset_check import check_inputs
[2]:
scenario = "ref"
suffix = "_02161113"
output_folder = f"{test_case_path}/output/{scenario}{suffix}"
solver = "highs" # ["gurobi", "xpress", "highs", "mosek"]
[3]:
config = read_config_file(file_path=test_case_path / "config.yaml")
model_parameters = build_input_parameters(config)
model_parameters = check_inputs(model_parameters)
WARNING:root:Missing variables in input dataset. Adding defaults:
storage_energy_to_power_ratio: nan (type=float64)
[4]:
model = build_model(model_parameters)
Flexibility demand sparsity: 25.00%
[5]:
model.solve(solver_name=solver)
save_solution(
model=model,
output_folder=output_folder,
model_parameters=model_parameters,
)
Running HiGHS 1.13.1 (git hash: n/a): Copyright (c) 2026 under MIT licence terms
LP linopy-problem-xt84_im2 has 6404 rows; 5950 cols; 19368 nonzeros
Coefficient ranges:
Matrix [1e-02, 9e+05]
Cost [1e-02, 1e+00]
Bound [0e+00, 0e+00]
RHS [1e+00, 1e+03]
Presolving model
3052 rows, 3222 cols, 10840 nonzeros 0s
2196 rows, 2390 cols, 9480 nonzeros 0s
Dependent equations search running on 624 equations with time limit of 1000.00s
Dependent equations search removed 0 rows and 0 nonzeros in 0.00s (limit = 1000.00s)
2068 rows, 2334 cols, 8024 nonzeros 0s
Presolve reductions: rows 2068(-4336); columns 2334(-3616); nonzeros 8024(-11344)
Solving the presolved LP
Using dual simplex solver
Iteration Objective Infeasibilities num(sum)
0 5.6017000644e-02 Pr: 280(231040) 0.0s
809 3.8396016000e+05 Pr: 0(0) 0.0s
Performed postsolve
Solving the original LP from the solution after postsolve
Model name : linopy-problem-xt84_im2
Model status : Optimal
Simplex iterations: 809
Objective value : 3.8396016000e+05
P-D objective error : 1.5159798119e-16
HiGHS run time : 0.04
[6]:
solution = model.solution
print(model.solution)
<xarray.Dataset> Size: 77kB
Dimensions: (area: 2, hour: 10, resource: 4,
year_op: 4, conversion_tech: 5,
year_dec: 6, year_inv: 5,
storage_tech: 3, link: 2,
transport_tech: 3,
combined_tech: 1, mode: 2,
hour_type: 5)
Coordinates: (12/13)
* area (area) <U6 48B 'area_1' 'area_2'
* hour (hour) int64 80B 0 1 2 3 4 5 6 7 8 9
* resource (resource) <U11 176B 'electricity...
* year_op (year_op) int64 32B 2020 ... 2050
* conversion_tech (conversion_tech) <U12 240B 'elec...
* year_dec (year_dec) int64 48B 2020 ... 2070
... ...
* storage_tech (storage_tech) <U13 156B 'battery...
* link (link) <U6 48B 'link_1' 'link_2'
* transport_tech (transport_tech) <U16 192B 'big_m...
* combined_tech (combined_tech) <U15 60B 'electri...
* mode (mode) <U8 64B 'electric' 'fossil'
* hour_type (hour_type) <U3 60B 'HCE' ... 'P'
Data variables: (12/47)
operation_load_shedding_power (area, hour, resource, year_op) float64 3kB ...
operation_spillage_power (area, hour, resource, year_op) float64 3kB ...
operation_load_shedding_costs (area, year_op) float64 64B -0.0 ...
operation_spillage_costs (area, year_op) float64 64B -0.0 ...
annualised_totex (area, year_op) float64 64B 4.831...
operation_conversion_power_capacity (area, conversion_tech, year_op) float64 320B ...
... ...
operation_turpe_contract_power (area, hour_type, year_op) float64 320B ...
operation_turpe_variable_costs (area, year_op) float64 64B 0.0 ....
operation_turpe_fixed_costs (area, year_op) float64 64B -0.0 ...
operation_flexibility_demand (area, resource, hour, year_op) float64 3kB ...
operation_flexibility_displaced_demand (area, resource, hour, year_op) float64 3kB ...
operation_flexibility_costs (area, resource, year_op) float64 256B ...
[7]:
dual = model.dual
print(dual)
<xarray.Dataset> Size: 97kB
Dimensions: (area: 2,
hour: 10,
resource: 4,
year_op: 4,
conversion_tech: 5,
year_inv: 5,
...
year_dec: 6,
link: 2,
transport_tech: 3,
combined_tech: 1,
hour_type: 5,
group: 5)
Coordinates: (12/13)
* area (area) <U6 48B ...
* hour (hour) int64 80B ...
* resource (resource) <U11 176B ...
* year_op (year_op) int64 32B ...
* conversion_tech (conversion_tech) <U12 240B ...
* year_inv (year_inv) int64 40B ...
... ...
* year_dec (year_dec) int64 48B ...
* link (link) <U6 48B ...
* transport_tech (transport_tech) <U16 192B ...
* combined_tech (combined_tech) <U15 60B ...
* hour_type (hour_type) <U3 60B ...
* group (group) int64 40B ...
Data variables: (12/95)
operation_adequacy_constraint (area, hour, resource, year_op) float64 3kB ...
operation_load_shedding_max_constraint (area, hour, resource, year_op) float64 3kB ...
operation_spillage_max_constraint (area, hour, resource, year_op) float64 3kB ...
operation_load_shedding_costs_def (area, year_op) float64 64B ...
operation_spillage_costs_def (area, year_op) float64 64B ...
annualised_totex_def (area, year_op) float64 64B ...
... ...
operation_flexibility_conservation_2h_hydrogen_area_2_2020 (group) float64 40B ...
operation_flexibility_conservation_2h_hydrogen_area_2_2030 (group) float64 40B ...
operation_flexibility_conservation_2h_hydrogen_area_2_2040 (group) float64 40B ...
operation_flexibility_conservation_2h_hydrogen_area_2_2050 (group) float64 40B ...
operation_flexibility_displaced_demand_def (area, resource, hour, year_op) float64 3kB ...
operation_flexibility_costs_def (area, resource, year_op, hour) float64 3kB ...
[8]:
prices = dual.operation_adequacy_constraint
df = prices.to_dataframe(name="value").reset_index()
print(df)
area hour resource year_op value
0 area_1 0 electricity 2020 -0.0
1 area_1 0 electricity 2030 -0.0
2 area_1 0 electricity 2040 -0.0
3 area_1 0 electricity 2050 -0.0
4 area_1 0 heat 2020 876.0
.. ... ... ... ... ...
315 area_2 9 hydrogen 2050 12.5
316 area_2 9 methane 2020 -0.0
317 area_2 9 methane 2030 -0.0
318 area_2 9 methane 2040 -0.0
319 area_2 9 methane 2050 -0.0
[320 rows x 5 columns]
[9]:
unique_resources = df["resource"].unique()
unique_years = df["year_op"].unique()
unique_areas = df["area"].unique()
colors = plt.cm.tab10(range(len(unique_areas)))
area_color_map = {area: colors[i] for i, area in enumerate(unique_areas)}
fig, axes = plt.subplots(
len(unique_resources),
len(unique_years),
figsize=(15, 10),
sharex=True,
sharey="row",
)
for i, resource in enumerate(unique_resources):
for j, year in enumerate(unique_years):
ax = axes[i, j]
subset = df[(df["resource"] == resource) & (df["year_op"] == year)]
for area in unique_areas:
area_data = subset[subset["area"] == area]
if not area_data.empty:
ax.plot(
area_data["hour"],
area_data["value"],
marker="o",
color=area_color_map[area],
)
if i == 0:
ax.set_title(f"Year: {year}")
if j == 0:
ax.set_ylabel(f"Resource: {resource}")
if i == len(unique_resources) - 1:
ax.set_xlabel("Hour")
legend_elements = [
plt.Line2D(
[0],
[0],
marker="o",
color=color,
linestyle="",
markersize=8,
label=area,
)
for area, color in area_color_map.items()
]
fig.legend(
handles=legend_elements,
loc="upper center",
ncol=len(unique_areas),
title="Area",
)
plt.tight_layout(rect=(0, 0, 1, 0.95))
plt.show()
[ ]: