Tutorial: OBZ case study

Tutorial for the OBZ case study as an example of the full workflow of Tulipa.

We are basing ourselves on the Tulipa data pipeline/workflow. To help us navigate this workflow, we'll reproduce the diagram from the link above here. For more details on the steps of the workflow, check the original link, or follow the tutorial.

Install packages

To follow this tutorial, you need to install some packages:

1. Open julia in the folder that you will be working
2. In the julia terminal, press ]. You should see pkg>
3. Activate you local environment using activate .
4. Install packages with add <Package1>,<Package2>, etc. You shouldn't need to specify the versions, if this tutorial is up-to-date.

These are the installed packages and their versions:

using Pkg
Pkg.status()

Status `~/work/TulipaEnergyModel.jl/TulipaEnergyModel.jl/docs/Project.toml`
  [336ed68f] CSV v0.10.15
  [a93c6f00] DataFrames v1.7.0
  [b4f34e82] Distances v0.10.12
  [e30172f5] Documenter v1.11.4
  [d2f5444f] DuckDB v1.2.2
  [60bf3e95] GLPK v1.2.1
  [87dc4568] HiGHS v1.17.0
  [682c06a0] JSON v0.21.4
  [4076af6c] JuMP v1.26.0
  [16fef848] LiveServer v1.5.0
  [bac558e1] OrderedCollections v1.8.1
  [91a5bcdd] Plots v1.40.13
  [314fac8b] TulipaClustering v0.4.1
  [5d7bd171] TulipaEnergyModel v0.16.0 `~/work/TulipaEnergyModel.jl/TulipaEnergyModel.jl`
  [7b3808b7] TulipaIO v0.5.0

External source

For this tutorial, we'll use the OBZ data. Download it from the Zenodo Link and store it in a folder.

Check https://github.com/TulipaEnergy/Tulipa-OBZ-CaseStudy for more information.

These are the files that we are working with:

# user_input_dir should point to the folder where the data was downloaded and extracted
readdir(user_input_dir)

20-element Vector{String}:
 "assets-consumer-basic-data.csv"
 "assets-consumer-yearly-data.csv"
 "assets-conversion-basic-data.csv"
 "assets-conversion-yearly-data.csv"
 "assets-hub-basic-data.csv"
 "assets-hub-yearly-data.csv"
 "assets-producer-basic-data.csv"
 "assets-producer-yearly-data.csv"
 "assets-profiles.csv"
 "assets-storage-basic-data.csv"
 "assets-storage-min-max-reservoir-level-profiles.csv"
 "assets-storage-yearly-data.csv"
 "flows-assets-connections-basic-data.csv"
 "flows-assets-connections-yearly-data.csv"
 "flows-transport-assets-basic-data.csv"
 "flows-transport-assets-yearly-data.csv"
 "inputs-description.toml"
 "min-max-reservoir-levels.csv"
 "profiles.csv"
 "year-data.csv"

For the Tulipa workflow, we will need to transform some of this data into a specific format. This can be done externally in whatever tools you are already comfortable with, or through Julia via DuckDB and Tulipa's convenience functions.

Create connection

Once we are done manipulating the data externally, it is time to create a DuckDB connection.

You can create a connection storing the DB locally, or keep everything in-memory only. Let's assume you want to store the DB, otherwise you can just remove the argument "obz.db".

using DuckDB: DBInterface, DuckDB

# We are staring from a fresh `obz.db` file
connection = DBInterface.connect(DuckDB.DB, "obz.db")

DuckDB.DB("obz.db")

We will be performing various queries with DuckDB. To format them nicely, we can wrap the results in a DataFrame:

using DataFrames: DataFrame

nice_query(str) = DataFrame(DuckDB.query(connection, str))

nice_query (generic function with 1 method)

Load data

Once we are done manipulating the data externally, it is time to load it into the DuckDB connection.

Note that this doesn't have to be the Tulipa-specific data. It can be whatever data you prefer to manipulate via Julia/DuckDB, instead of externally.

We can load them manually with DuckDB, but we also have a convenience function:

using TulipaIO: TulipaIO

TulipaIO.read_csv_folder(
    connection,
    user_input_dir,
    replace_if_exists = true,
)

# The first 5 tables
nice_query("SELECT table_name FROM duckdb_tables() LIMIT 5")

Data processing for instance data with DuckDB/TulipaIO

As we mentioned before, you can process your data externally and then load it. But you can also use Julia and DuckDB to process the data.

This step is required to prepare the data for TulipaClustering for the clustering of the profile data.

We need a single profiles table with 4 columns:

• profile_name
• year
• timestep
• value

Instead, we have the profiles data in the profiles table, which looks something like the following, but with many more columns:

nice_query("SELECT year, timestep, * LIKE 'NL_%' FROM profiles LIMIT 5")

The total number of columns in the profiles table:

nice_query("SELECT COUNT(*) FROM duckdb_columns() WHERE table_name = 'profiles'")

Notice that these are all hourly profiles for the whole year:

nice_query("SELECT year, MAX(timestep) FROM profiles GROUP BY year")

So we will transform both this table to long format:

using TulipaClustering: TulipaClustering

TulipaClustering.transform_wide_to_long!(connection, "profiles", "pivot_profiles")

DuckDB.query(
    connection,
    "CREATE OR REPLACE TABLE profiles AS
    FROM pivot_profiles
    ORDER BY profile_name, year, timestep
    "
)

nice_query("SELECT COUNT(*) FROM profiles")

Just to showcase this, let's plot all NL_* profiles in the first 72 hours of the year:

using Plots

subtable = DuckDB.query(
    connection,
    "SELECT
        timestep,
        value,
        profile_name,
    FROM profiles
    WHERE
        profile_name LIKE 'NL_%'
        AND year=2050
        AND timestep <= 72 -- Just 72 hours
    ORDER BY timestep
    ",
)
df = DataFrame(subtable)
plot(df.timestep, df.value, group=df.profile_name)

Cluster into representative periods using TulipaClustering

Instead of working with the full time horizon of 8760 hours, we will cluster the profiles using TulipaClustering.

You can tweak around the number of representative periods and period duration (and naturally other parameters). This is the configuration that we'll use:

using Distances: SqEuclidean

## Data for clustering
clustering_params = (
    num_rep_periods = 3,    # number of representative periods
    period_duration = 24,   # hours of the representative period
    method = :k_means,
    distance = SqEuclidean(),
    ## Data for weight fitting
    weight_type = :convex,
    tol = 1e-2,
)

(num_rep_periods = 3, period_duration = 24, method = :k_means, distance = Distances.SqEuclidean(0.0), weight_type = :convex, tol = 0.01)

PS. We chose to define our clustering_params as a NamedTuple in Julia, but that is completely optional and you can use whatever structure suits your case.

using Random: Random
Random.seed!(123)
TulipaClustering.cluster!(
    connection,
    clustering_params.period_duration,  # Required
    clustering_params.num_rep_periods;  # Required
    clustering_params.method,           # Optional
    clustering_params.distance,         # Optional
    clustering_params.weight_type,      # Optional
    clustering_params.tol,              # Optional
);

TulipaClustering.ClusteringResult(11736×5 DataFrame
   Row │ rep_period  timestep  year   profile_name     value
       │ Int64       Int64     Int64  String           Float64
───────┼────────────────────────────────────────────────────────
     1 │          1         1   2050  AT_E_Demand      0.433742
     2 │          1         2   2050  AT_E_Demand      0.413336
     3 │          1         3   2050  AT_E_Demand      0.400568
     4 │          1         4   2050  AT_E_Demand      0.404315
     5 │          1         5   2050  AT_E_Demand      0.447875
     6 │          1         6   2050  AT_E_Demand      0.490754
     7 │          1         7   2050  AT_E_Demand      0.555144
     8 │          1         8   2050  AT_E_Demand      0.592532
   ⋮   │     ⋮          ⋮        ⋮           ⋮            ⋮
 11730 │          3        18   2050  UK_Wind_Onshore  0.770287
 11731 │          3        19   2050  UK_Wind_Onshore  0.769074
 11732 │          3        20   2050  UK_Wind_Onshore  0.769339
 11733 │          3        21   2050  UK_Wind_Onshore  0.769764
 11734 │          3        22   2050  UK_Wind_Onshore  0.768185
 11735 │          3        23   2050  UK_Wind_Onshore  0.766419
 11736 │          3        24   2050  UK_Wind_Onshore  0.75759
                                              11721 rows omitted, sparse([46, 47, 48, 49, 51, 53, 62, 67, 69, 70  …  356, 357, 358, 359, 360, 361, 362, 363, 364, 365], [1, 1, 1, 1, 1, 1, 1, 1, 1, 1  …  3, 3, 3, 3, 3, 3, 3, 3, 3, 3], [0.14481082337653622, 0.13121132447100464, 0.04959948680279389, 0.05227205608976673, 0.11217960913380766, 0.03744057362266463, 0.014585541277927053, 0.04976782552689989, 0.47021910388116805, 0.09138795402559592  …  1.0, 1.0, 1.0, 1.0, 0.9894463485115309, 0.8167218510506219, 0.7836474990520129, 0.8652566754986826, 1.0, 1.0], 365, 3), [0.561845156 0.715936085 … 0.722342321 0.692761049; 0.533957227 0.690376332 … 0.69299084 0.666338922; … ; 0.944921215 0.84162672 … 0.572282407 0.888713937; 0.943940536 0.870985063 … 0.602038058 0.870050813], [0.433741974188679 0.4888906324253732 0.6131286770479999; 0.4133362753867923 0.4679652152985074 0.5866607760399997; … ; 0.5037757928113207 0.5261235292985074 0.7664189462400003; 0.5021113885377356 0.5181840313432834 0.757589973752], TulipaClustering.AuxiliaryClusteringData([:timestep, :year, :profile_name], 24, 24, 365, nothing))

These are the tables created by TulipaClustering:

• rep_periods_data
• rep_periods_mapping
• profiles_rep_periods
• timeframe_data

Prepare data for TulipaEnergyModel's format

Now the fun part starts. We need to create specific tables for Tulipa using our current tables. Again, we remind you that you can create most of these files externally, i.e., you don't have to use DuckDB to join them here. However, defining the workflow in a programmatic way makes it easier to reproduce it in the future.

We have to define a minimum set of columns for each table, and then the remaining columns will be filled with defaults. Some columns cannot contain missing values (such as the asset or year columns in most tables). For other columns, missing values will be filled with the columns' default.

Year data

This data is already correct in the case study and contains a single year.

nice_query("FROM year_data")

Assets

First, let's join all assets' basic data. This table goes directly into the asset table for Tulipa.

DuckDB.query(
    connection,
    "CREATE TABLE asset AS
    SELECT
        name AS asset,
        type,
        capacity,
        capacity_storage_energy,
        is_seasonal,
    FROM (
        FROM assets_consumer_basic_data
        UNION BY NAME
        FROM assets_conversion_basic_data
        UNION BY NAME
        FROM assets_hub_basic_data
        UNION BY NAME
        FROM assets_producer_basic_data
        UNION BY NAME
        FROM assets_storage_basic_data
    )
    ORDER BY asset
    ",
)

nice_query("FROM asset ORDER BY random() LIMIT 5")

Similarly, we join the assets' yearly data:

DuckDB.query(
    connection,
    "CREATE TABLE t_asset_yearly AS
    FROM (
        FROM assets_consumer_yearly_data
        UNION BY NAME
        FROM assets_conversion_yearly_data
        UNION BY NAME
        FROM assets_hub_yearly_data
        UNION BY NAME
        FROM assets_producer_yearly_data
        UNION BY NAME
        FROM assets_storage_yearly_data
    )
    ",
)

nice_query("FROM t_asset_yearly ORDER BY random() LIMIT 5")

Then, the t_asset_yearly table is used to create the three other asset tables that Tulipa requires:

DuckDB.query(
    connection,
    "CREATE TABLE asset_commission AS
    SELECT
        name AS asset,
        year AS commission_year,
    FROM t_asset_yearly
    ORDER by asset
    "
)

DuckDB.query(
    connection,
    "CREATE TABLE asset_milestone AS
    SELECT
        name AS asset,
        year AS milestone_year,
        peak_demand,
        initial_storage_level,
        storage_inflows,
    FROM t_asset_yearly
    ORDER by asset
    "
)

DuckDB.query(
    connection,
    "CREATE TABLE asset_both AS
    SELECT
        name AS asset,
        year AS milestone_year,
        year AS commission_year, -- Yes, it is the same year twice with different names because it's not a multi-year problem
        initial_units,
        initial_storage_units,
    FROM t_asset_yearly
    ORDER by asset
    "
)

(Count = [391],)

Here is an example of one of these tables:

nice_query("FROM asset_both WHERE initial_storage_units > 0 LIMIT 5")

Flows

We repeat the steps above for flows:

DuckDB.query(
    connection,
    "CREATE TABLE flow AS
    SELECT
        from_asset,
        to_asset,
        carrier,
        capacity,
        is_transport,
    FROM (
        FROM flows_assets_connections_basic_data
        UNION BY NAME
        FROM flows_transport_assets_basic_data
    )
    ORDER BY from_asset, to_asset
    ",
)

DuckDB.query(
    connection,
    "CREATE TABLE t_flow_yearly AS
    FROM (
        FROM flows_assets_connections_yearly_data
        UNION BY NAME
        FROM flows_transport_assets_yearly_data
    )
    ",
)

DuckDB.query(
    connection,
    "CREATE TABLE flow_commission AS
    SELECT
        from_asset,
        to_asset,
        year AS commission_year,
        efficiency,
    FROM t_flow_yearly
    ORDER by from_asset, to_asset
    "
)

DuckDB.query(
    connection,
    "CREATE TABLE flow_milestone AS
    SELECT
        from_asset,
        to_asset,
        year AS milestone_year,
        variable_cost,
    FROM t_flow_yearly
    ORDER by from_asset, to_asset
    "
)

DuckDB.query(
    connection,
    "CREATE TABLE flow_both AS
    SELECT
        from_asset,
        to_asset,
        year AS milestone_year,
        year AS commission_year,
        initial_export_units,
        initial_import_units,
    FROM t_flow_yearly
    ORDER by from_asset, to_asset
    "
)

(Count = [501],)

Assets profiles

The assets_profiles table already exists, so we only need to create assets_timeframe_profiles. Since all the data is already in assets_storage_min_max_reservoir_level_profiles, we just copy it over.

DuckDB.query(
    connection,
    "CREATE TABLE assets_timeframe_profiles AS
    FROM assets_storage_min_max_reservoir_level_profiles
    ORDER BY asset, commission_year, profile_name
    ",
)

(Count = [12],)

Partitions

The OBZ table uses only uniform time partitions, which makes it easy to create the necessary tables.

For the assets_rep_periods_partitions, we simply have to copy the partition given by the assets' yearly data for each representative period given rep_periods_data.rep_period and attach a specification = 'uniform' to that table.

DuckDB.query(
    connection,
    "CREATE TABLE assets_rep_periods_partitions AS
    SELECT
        t.name AS asset,
        t.year,
        t.partition AS partition,
        rep_periods_data.rep_period,
        'uniform' AS specification,
    FROM t_asset_yearly AS t
    LEFT JOIN rep_periods_data
        ON t.year = rep_periods_data.year
    ORDER BY asset, t.year, rep_period
    ",
)

(Count = [1173],)

For the flows_rep_periods_partitions, we need to also compute the expected partition value, which will follow a simple formula. Given a flow (from_asset, to_asset), we look at the partition of both from_asset and to_asset. If the flow is a transport flow, we use the maximum between the partitions of from_asset and to_asset. Otherwise, we use the minimum between these two.

DuckDB.query(
    connection,
    "CREATE TABLE flows_rep_periods_partitions AS
    SELECT
        flow.from_asset,
        flow.to_asset,
        t_from.year,
        t_from.rep_period,
        'uniform' AS specification,
        IF(
            flow.is_transport,
            greatest(t_from.partition::int, t_to.partition::int),
            least(t_from.partition::int, t_to.partition::int)
        ) AS partition,
    FROM flow
    LEFT JOIN assets_rep_periods_partitions AS t_from
        ON flow.from_asset = t_from.asset
    LEFT JOIN assets_rep_periods_partitions AS t_to
        ON flow.to_asset = t_to.asset
        AND t_from.year = t_to.year
        AND t_from.rep_period = t_to.rep_period
    ",
)

(Count = [1503],)

Timeframe profiles

For the timeframe profiles, we'll use the other profiles table that we haven't touched yet: min_max_reservoir_levels. As with the other profiles, we will first pivot that table to have it in long format. However, we do not cluster these profiles, since the representative periods are already computed. Instead, we will create a temporary table (cte_split_profiles) that converts the timestep that goes from 1 to 8760 into two columns: period, from to 1 to 365 (days) and timestep, from 1 to 24 (hours).

Finally, the timeframe profiles are computed with the average over period, i.e., each value of a given timeframe profile in a period is the average of 24 hours of the original profile.

TulipaClustering.transform_wide_to_long!(
    connection,
    "min_max_reservoir_levels",
    "pivot_min_max_reservoir_levels",
)

period_duration = clustering_params.period_duration

DuckDB.query(
    connection,
    "
    CREATE TABLE profiles_timeframe AS
    WITH cte_split_profiles AS (
        SELECT
            profile_name,
            year,
            1 + (timestep - 1) // $period_duration  AS period,
            1 + (timestep - 1)  % $period_duration AS timestep,
            value,
        FROM pivot_min_max_reservoir_levels
    )
    SELECT
        cte_split_profiles.profile_name,
        cte_split_profiles.year,
        cte_split_profiles.period,
        AVG(cte_split_profiles.value) AS value, -- Computing the average aggregation
    FROM cte_split_profiles
    GROUP BY
        cte_split_profiles.profile_name,
        cte_split_profiles.year,
        cte_split_profiles.period
    ORDER BY
        cte_split_profiles.profile_name,
        cte_split_profiles.year,
        cte_split_profiles.period
    ",
)

(Count = [4380],)

Populate with defaults

Finally, in many cases, you will need to complete the missing columns with additional information. To simplify this process, we created the populate_with_defaults! function. Please read TulipaEnergyModel's populate with default section for a complete picture.

Here is the before of one of the tables:

nice_query("FROM asset_both LIMIT 5")

using TulipaEnergyModel: TulipaEnergyModel as TEM

TEM.populate_with_defaults!(connection)

nice_query("FROM asset_both LIMIT 5")

Create internal tables for the model indices

energy_problem = TEM.EnergyProblem(connection)

EnergyProblem:
  - Model not created!
  - Model not solved!

Purely out of curiosity, here is the total number of tables that we have:

nice_query("SELECT COUNT(*) as num_tables, FROM duckdb_tables()")

Create model

Finally, we get to actually use the model.

optimizer_parameters = Dict(
    "output_flag" => true,
    "mip_rel_gap" => 0.0,
    "mip_feasibility_tolerance" => 1e-5,
)
TEM.create_model!(energy_problem; model_file_name, optimizer_parameters)

EnergyProblem:
  - Model created!
    - Number of variables: 34075
    - Number of constraints for variable bounds: 32221
    - Number of structural constraints: 72855
  - Model not solved!

Solve model

Last, but not least important, we solve the model:

TEM.solve_model!(energy_problem)

Running HiGHS 1.10.0 (git hash: fd8665394e): Copyright (c) 2025 HiGHS under MIT licence terms
LP   has 72855 rows; 34075 cols; 586771 nonzeros
Coefficient ranges:
  Matrix [9e-03, 6e+00]
  Cost   [1e+03, 5e+05]
  Bound  [4e-01, 5e+04]
  RHS    [1e-06, 9e+04]
Presolving model
18937 rows, 28752 cols, 524963 nonzeros  0s
Dependent equations search running on 18139 equations with time limit of 1000.00s
Dependent equations search removed 242 rows and 6492 nonzeros in 0.07s (limit = 1000.00s)
17876 rows, 25581 cols, 499670 nonzeros  0s
Presolve : Reductions: rows 17876(-54979); columns 25581(-8494); elements 499670(-87101)
Solving the presolved LP
Using EKK dual simplex solver - serial
  Iteration        Objective     Infeasibilities num(sum)
          0     0.0000000000e+00 Ph1: 0(0) 1s
      20603     9.3979319318e+06 Pr: 0(0); Du: 0(4.35755e-09) 2s
Solving the original LP from the solution after postsolve
Model status        : Optimal
Simplex   iterations: 20603
Objective value     :  9.3979319318e+06
Relative P-D gap    :  7.1351043888e-15
HiGHS run time      :          2.19

Store primal and dual solution

The primal and dual solutions are computed when saving the solution with save_solution!, as long as we don't change the default value of compute_duals. Here it is, explicitly:

TEM.save_solution!(energy_problem; compute_duals = true)

Now every variable indices table has a column solution, and every constraint has additional columns dual_*, depending on the constraints name.

Examples checking the primal and dual solutions

Select all variables of storage level at representative periods with value greater than 0 at the solution (show only first 5):

nice_query("SELECT *
    FROM var_storage_level_rep_period
    WHERE solution > 0
    LIMIT 5
")

Select all indices related to the balance storage at representative periods when both min_storage_level_rep_period_limit and max_storage_level_rep_period_limit have duals equal to 0.

nice_query("SELECT *
    FROM cons_balance_storage_rep_period
    WHERE dual_max_storage_level_rep_period_limit = 0
        AND dual_min_storage_level_rep_period_limit = 0
    LIMIT 5
")

Data processing for plots and dashboard

This part of the workflow is open for you to do whatever you need. In principle, you can skip this step and go straight to exporting the solution, and then perform your analysis of the solution outside of the DuckDB/Julia environment.

Here is an example of data processing using DuckDB and Julia.

The table

nice_query("
CREATE TEMP TABLE analysis_inter_storage_levels AS
SELECT
    var.id,
    var.asset,
    var.period_block_start as period,
    asset.capacity_storage_energy,
    var.solution / (
        IF(asset.capacity_storage_energy > 0, asset.capacity_storage_energy, 1)
    ) AS SoC,
FROM var_storage_level_over_clustered_year AS var
LEFT JOIN asset
    ON var.asset = asset.asset
")
nice_query("FROM analysis_inter_storage_levels LIMIT 5")

Create plots

Now, using the analysis tables from above, we can create plots in Julia:

using Plots

p = plot()
assets = ["ES_Hydro_Reservoir", "NO_Hydro_Reservoir", "FR_Hydro_Reservoir"]

df = nice_query("SELECT asset, period, SoC
    FROM analysis_inter_storage_levels
    WHERE asset in ('ES_Hydro_Reservoir', 'NO_Hydro_Reservoir', 'FR_Hydro_Reservoir')
")

plot!(
    df.period,          # x-axis
    df.SoC,             # y-axis
    group = df.asset,   # each asset is a different plot
    xlabel = "Period",
    ylabel = "Storage level [p.u.]",
    linewidth = 3,
    dpi = 600,
)

Export solution

Finally, we can export the solution to CSV files using the convenience function below:

mkdir("obz-outputs")
TEM.export_solution_to_csv_files("obz-outputs", energy_problem)
readdir("obz-outputs")

11-element Vector{String}:
 "cons_balance_consumer.csv"
 "cons_balance_conversion.csv"
 "cons_balance_hub.csv"
 "cons_balance_storage_over_clustered_year.csv"
 "cons_balance_storage_rep_period.csv"
 "cons_capacity_incoming_simple_method.csv"
 "cons_capacity_outgoing_simple_method.csv"
 "cons_transport_flow_limit_simple_method.csv"
 "var_flow.csv"
 "var_storage_level_over_clustered_year.csv"
 "var_storage_level_rep_period.csv"

Using DuckDB directly it is also possible to export to other formats, such as Parquet.