Tutorials

Here are some tutorials on how to use Tulipa.

• Tutorials
• • Basic example
  • • Run scenario
    • Manually running each step
    • Manually creating all structures without EnergyProblem
    • Change optimizer and specify parameters
    • Using the graph structure
  • Manipulating the solution
  • • The solution returned by solve_model
    • The solution inside the graph
    • The solution inside the dataframes object
    • Values of constraints and expressions
    • Writing the output to CSV
    • Plotting
  • Hydrothermal Dispatch example

Basic example

For our first example, let's use a tiny existing dataset. Inside the code for this package, you can find the folder test/inputs/Tiny, which includes all the files necessary to create a model and solve it.

The files inside the "Tiny" folder define the assets and flows data, their profiles, and their time resolution, as well as define the representative periods and which periods in the full problem formulation they represent.¹

For more details about these files, see Input.

¹ Ignore bad-assets-data.csv, which is used for testing.

Run scenario

To read all data from the Tiny folder, perform all necessary steps to create a model, and solve the model, run the following in a Julia terminal:

using TulipaEnergyModel

# input_dir should be the path to Tiny as a string (something like "test/inputs/Tiny")
energy_problem = run_scenario(input_dir)

EnergyProblem:
  - Time for reading the data (in seconds): 0.097166277
  - Model created!
    - Time for creating the model (in seconds): 0.296669415
    - Number of variables: 364
    - Number of constraints for variable bounds: 364
    - Number of structural constraints: 432
  - Model solved! 
    - Time for solving the model (in seconds): 0.710980233
    - Termination status: OPTIMAL
    - Objective value: 269238.4382375954

The energy_problem variable is of type EnergyProblem. For more details, see the documentation for that type or the section Structures.

That's all it takes to run a scenario! To learn about the data required to run your own scenario, see the Input section of How to Use.

Manually running each step

If we need more control, we can create the energy problem first, then the optimization model inside it, and finally ask for it to be solved.

using TulipaEnergyModel

# input_dir should be the path to Tiny
energy_problem = create_energy_problem_from_csv_folder(input_dir)

EnergyProblem:
  - Time for reading the data (in seconds): NaN
  - Model not created!
  - Model not solved!

The energy problem does not have a model yet:

energy_problem.model === nothing

true

To create the internal model, we call the function create_model!.

create_model!(energy_problem)
energy_problem.model

A JuMP Model
Minimization problem with:
Variables: 364
Objective function type: JuMP.AffExpr
`JuMP.AffExpr`-in-`MathOptInterface.EqualTo{Float64}`: 72 constraints
`JuMP.AffExpr`-in-`MathOptInterface.LessThan{Float64}`: 360 constraints
`JuMP.VariableRef`-in-`MathOptInterface.GreaterThan{Float64}`: 364 constraints
`JuMP.VariableRef`-in-`MathOptInterface.LessThan{Float64}`: 1 constraint
`JuMP.VariableRef`-in-`MathOptInterface.Integer`: 4 constraints
Model mode: AUTOMATIC
CachingOptimizer state: NO_OPTIMIZER
Solver name: No optimizer attached.
Names registered in the model: assets_investment, assets_profile_times_capacity_in, assets_profile_times_capacity_out, consumer_balance, conversion_balance, energy_limit, flow, flows_investment, hub_balance, incoming_flow_highest_in_out_resolution, incoming_flow_highest_in_resolution, incoming_flow_lowest_resolution, incoming_flow_lowest_storage_resolution_intra_rp, incoming_flow_storage_inter_rp_balance, max_input_flows_limit, max_output_flows_limit, max_storage_level_inter_rp_limit, max_storage_level_intra_rp_limit, max_transport_flow_limit, min_storage_level_inter_rp_limit, min_storage_level_intra_rp_limit, min_transport_flow_limit, outgoing_flow_highest_in_out_resolution, outgoing_flow_highest_out_resolution, outgoing_flow_lowest_resolution, outgoing_flow_lowest_storage_resolution_intra_rp, outgoing_flow_storage_inter_rp_balance, storage_level_inter_rp, storage_level_intra_rp

The model has not been solved yet, which can be verified through the solved flag inside the energy problem:

energy_problem.solved

false

Finally, we can solve the model:

solution = solve_model!(energy_problem)

TulipaEnergyModel.Solution(Dict(:ocgt => 1.1368683772161603e-15, :wind => 35.0, :ccgt => 2.0, :solar => 44.99999994237288), Dict{Tuple{Symbol, Symbol}, Float64}(), Float64[], Float64[], [0.0, 1.1368683772161603e-13, 1.1368683772161603e-12, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], 269238.4382388837, Dict(:hub_balance => [], :consumer_balance => [3.6500000000000004, 3.6500000000000004, 3.6500000000000004, 3.6500000000000004, 3.6500000000000004, 3.6500000000000004, 3.6500000000000004, 3.6500000000000004, 3.6500000000000004, 3.6500000000000004  …  9.125, 9.125, 9.125, 9.125, 9.125, 9.125, 9.125, 9.125, 9.125, 9.125]))

The solution is included in the individual assets and flows, but for completeness, we return the full solution object, also defined in the Structures section.

In particular, the objective value and the termination status are also included in the energy problem:

energy_problem.objective_value, energy_problem.termination_status

(269238.4382375954, MathOptInterface.OPTIMAL)

Manually creating all structures without EnergyProblem

For additional control, it might be desirable to use the internal structures of EnergyProblem directly. This can be error-prone, but it is slightly more efficient. The full description for these structures can be found in Structures.

using TulipaEnergyModel

# input_dir should be the path to Tiny
graph, representative_periods, timeframe = create_graph_and_representative_periods_from_csv_folder(input_dir)

(Meta graph based on a Graphs.SimpleGraphs.SimpleDiGraph{Int64} with vertex labels of type Symbol, vertex metadata of type GraphAssetData, edge metadata of type GraphFlowData, graph metadata given by nothing, and default weight nothing, RepresentativePeriod[RepresentativePeriod(Dict(352 => 1.0, 312 => 1.0, 81 => 1.0, 341 => 1.0, 141 => 1.0, 32 => 1.0, 291 => 1.0, 31 => 1.0, 282 => 1.0, 272 => 1.0…), 73.0, 1:24, 1.0), RepresentativePeriod(Dict(56 => 1.0, 35 => 1.0, 6 => 1.0, 67 => 1.0, 234 => 1.0, 215 => 1.0, 73 => 1.0, 164 => 1.0, 115 => 1.0, 153 => 1.0…), 182.5, 1:24, 1.0), RepresentativePeriod(Dict(110 => 1.0, 60 => 1.0, 220 => 1.0, 30 => 1.0, 268 => 1.0, 308 => 1.0, 219 => 1.0, 319 => 1.0, 320 => 1.0, 348 => 1.0…), 109.5, 1:24, 1.0)], Timeframe(365, 366×3 DataFrame
 Row │ period  rep_period  weight
     │ Int64   Int64       Float64
─────┼─────────────────────────────
   1 │      1           1      1.0
   2 │      2           1      1.0
   3 │      3           2      1.0
   4 │      4           2      1.0
   5 │      5           2      1.0
   6 │      6           2      1.0
   7 │      7           2      1.0
   8 │      8           3      1.0
  ⋮  │   ⋮         ⋮          ⋮
 360 │    360           3      1.0
 361 │    361           1      1.0
 362 │    362           2      1.0
 363 │    363           2      1.0
 364 │    364           3      1.0
 365 │    365           2      0.5
 366 │    365           3      0.5
                   351 rows omitted))

We also need a time partition for the constraints to create the model. Creating an energy problem automatically computes this data, but since we are doing it manually, we need to calculate it ourselves.

constraints_partitions = compute_constraints_partitions(graph, representative_periods)

Dict{Symbol, Dict{Tuple{Symbol, Int64}, Vector{UnitRange{Int64}}}} with 5 entries:
  :lowest                        => Dict((:ccgt, 3)=>[1:1, 2:2, 3:3, 4:4, 5:5, …
  :highest_in                    => Dict()
  :lowest_storage_level_intra_rp => Dict()
  :highest_in_out                => Dict((:demand, 1)=>[1:1, 2:2, 3:3, 4:4, 5:5…
  :highest_out                   => Dict((:ccgt, 3)=>[1:1, 2:2, 3:3, 4:4, 5:5, …

The constraints_partitions has two dictionaries with the keys :lowest_resolution and :highest_resolution. The lowest resolution dictionary is mainly used to create the constraints for energy balance, whereas the highest resolution dictionary is mainly used to create the capacity constraints in the model.

Finally, we also need dataframes that store the linearized indexes of the variables.

dataframes = construct_dataframes(graph, representative_periods, constraints_partitions, timeframe)

Dict{Symbol, DataFrames.DataFrame} with 7 entries:
  :lowest                        => 360×4 DataFrame…
  :highest_in                    => 0×4 DataFrame…
  :lowest_storage_level_intra_rp => 0×4 DataFrame…
  :highest_in_out                => 72×4 DataFrame…
  :highest_out                   => 360×4 DataFrame…
  :flows                         => 360×6 DataFrame…
  :storage_level_inter_rp        => 0×3 DataFrame…

Now we can compute the model.

model = create_model(graph, representative_periods, dataframes, timeframe)

A JuMP Model
Minimization problem with:
Variables: 364
Objective function type: JuMP.AffExpr
`JuMP.AffExpr`-in-`MathOptInterface.EqualTo{Float64}`: 72 constraints
`JuMP.AffExpr`-in-`MathOptInterface.LessThan{Float64}`: 360 constraints
`JuMP.VariableRef`-in-`MathOptInterface.GreaterThan{Float64}`: 364 constraints
`JuMP.VariableRef`-in-`MathOptInterface.LessThan{Float64}`: 1 constraint
`JuMP.VariableRef`-in-`MathOptInterface.Integer`: 4 constraints
Model mode: AUTOMATIC
CachingOptimizer state: NO_OPTIMIZER
Solver name: No optimizer attached.
Names registered in the model: assets_investment, assets_profile_times_capacity_in, assets_profile_times_capacity_out, consumer_balance, conversion_balance, energy_limit, flow, flows_investment, hub_balance, incoming_flow_highest_in_out_resolution, incoming_flow_highest_in_resolution, incoming_flow_lowest_resolution, incoming_flow_lowest_storage_resolution_intra_rp, incoming_flow_storage_inter_rp_balance, max_input_flows_limit, max_output_flows_limit, max_storage_level_inter_rp_limit, max_storage_level_intra_rp_limit, max_transport_flow_limit, min_storage_level_inter_rp_limit, min_storage_level_intra_rp_limit, min_transport_flow_limit, outgoing_flow_highest_in_out_resolution, outgoing_flow_highest_out_resolution, outgoing_flow_lowest_resolution, outgoing_flow_lowest_storage_resolution_intra_rp, outgoing_flow_storage_inter_rp_balance, storage_level_inter_rp, storage_level_intra_rp

Finally, we can compute the solution.

solution = solve_model(model)

TulipaEnergyModel.Solution(Dict(:ocgt => 1.1368683772161603e-15, :wind => 35.0, :ccgt => 2.0, :solar => 44.99999994237288), Dict{Tuple{Symbol, Symbol}, Float64}(), Float64[], Float64[], [0.0, 1.1368683772161603e-13, 1.1368683772161603e-12, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], 269238.4382388837, Dict(:hub_balance => [], :consumer_balance => [3.6500000000000004, 3.6500000000000004, 3.6500000000000004, 3.6500000000000004, 3.6500000000000004, 3.6500000000000004, 3.6500000000000004, 3.6500000000000004, 3.6500000000000004, 3.6500000000000004  …  9.125, 9.125, 9.125, 9.125, 9.125, 9.125, 9.125, 9.125, 9.125, 9.125]))

or, if we want to store the flow, storage_level_intra_rp, and storage_level_inter_rp optimal value in the dataframes:

solution = solve_model!(dataframes, model)

TulipaEnergyModel.Solution(Dict(:ocgt => 1.1368683772161603e-15, :wind => 35.0, :ccgt => 2.0, :solar => 44.99999994237288), Dict{Tuple{Symbol, Symbol}, Float64}(), Float64[], Float64[], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], 269238.4382375954, Dict(:hub_balance => [], :consumer_balance => [3.6500000000000004, 3.6500000000000004, 3.6500000000000004, 3.6500000000000004, 3.6500000000000004, 3.6500000000000004, 3.6500000000000004, 3.6500000000000004, 3.6500000000000004, 3.6500000000000004  …  9.125, 9.125, 9.125, 9.125, 9.125, 9.125, 9.125, 9.125, 9.125, 9.125]))

This solution structure is the same as the one returned when using an EnergyProblem.

Change optimizer and specify parameters

By default, the model is solved using the HiGHS optimizer (or solver). To change this, we can give the functions run_scenario, solve_model, or solve_model! a different optimizer.

For instance, we run the GLPK optimizer below:

using TulipaEnergyModel, GLPK

energy_problem = run_scenario(input_dir, optimizer = GLPK.Optimizer)

EnergyProblem:
  - Time for reading the data (in seconds): 0.003046825
  - Model created!
    - Time for creating the model (in seconds): 0.005618032
    - Number of variables: 364
    - Number of constraints for variable bounds: 364
    - Number of structural constraints: 432
  - Model solved! 
    - Time for solving the model (in seconds): 1.762546929
    - Termination status: OPTIMAL
    - Objective value: 269238.4382417078

or

using GLPK

solution = solve_model!(energy_problem, GLPK.Optimizer)

TulipaEnergyModel.Solution(Dict(:ocgt => 1.1368683772161603e-15, :wind => 35.0, :ccgt => 2.0, :solar => 44.99999994237288), Dict{Tuple{Symbol, Symbol}, Float64}(), Float64[], Float64[], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], 269238.4382375952, Dict(:hub_balance => [], :consumer_balance => [3.6500000000000004, 3.6500000000000004, 3.6500000000000004, 3.6500000000000004, 3.6500000000000004, 3.6500000000000004, 3.6500000000000004, 3.6500000000000004, 3.6500000000000004, 3.6500000000000004  …  9.125, 9.125, 9.125, 9.125, 9.125, 9.125, 9.125, 9.125, 9.125, 9.125]))

or

using GLPK

solution = solve_model(model, GLPK.Optimizer)

TulipaEnergyModel.Solution(Dict(:ocgt => 1.1368683772161603e-15, :wind => 35.0, :ccgt => 2.0, :solar => 44.99999994237288), Dict{Tuple{Symbol, Symbol}, Float64}(), Float64[], Float64[], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], 269238.4382375952, Dict(:hub_balance => [], :consumer_balance => [3.6500000000000004, 3.6500000000000004, 3.6500000000000004, 3.6500000000000004, 3.6500000000000004, 3.6500000000000004, 3.6500000000000004, 3.6500000000000004, 3.6500000000000004, 3.6500000000000004  …  9.125, 9.125, 9.125, 9.125, 9.125, 9.125, 9.125, 9.125, 9.125, 9.125]))

Notice that, in any of these cases, we need to explicitly add the GLPK package ourselves and add using GLPK before using GLPK.Optimizer.

In any of these cases, default parameters for the GLPK optimizer are used, which you can query using default_parameters. You can pass a dictionary using the keyword argument parameters to change the defaults. For instance, in the example below, we change the maximum allowed runtime for GLPK to be 1 seconds, which will most likely cause it to fail to converge in time.

using TulipaEnergyModel, GLPK

parameters = Dict("tm_lim" => 1)
energy_problem = run_scenario(input_dir, optimizer = GLPK.Optimizer, parameters = parameters)
energy_problem.termination_status

TIME_LIMIT::TerminationStatusCode = 12

For the complete list of parameters, check your chosen optimizer.

These parameters can also be passed via a file. See the read_parameters_from_file function for more details.

Using the graph structure

Read about the graph structure in the Graph section first.

We will use the graph created above for the "Tiny" dataset.

The first thing that we can do is access all assets. They are the labels of the graph and can be accessed via the MetaGraphsNext API:

using MetaGraphsNext
# Accessing assets
labels(graph)

Base.Generator{Base.OneTo{Int64}, MetaGraphsNext.var"#20#21"{MetaGraphsNext.MetaGraph{Int64, Graphs.SimpleGraphs.SimpleDiGraph{Int64}, Symbol, GraphAssetData, GraphFlowData, Nothing, Nothing, Nothing}}}(MetaGraphsNext.var"#20#21"{MetaGraphsNext.MetaGraph{Int64, Graphs.SimpleGraphs.SimpleDiGraph{Int64}, Symbol, GraphAssetData, GraphFlowData, Nothing, Nothing, Nothing}}(Meta graph based on a Graphs.SimpleGraphs.SimpleDiGraph{Int64} with vertex labels of type Symbol, vertex metadata of type GraphAssetData, edge metadata of type GraphFlowData, graph metadata given by nothing, and default weight nothing), Base.OneTo(6))

Notice that the result is a generator, so if we want the actual results, we have to collect it:

labels(graph) |> collect

6-element Vector{Symbol}:
 :ocgt
 :ccgt
 :wind
 :solar
 :ens
 :demand

To access the asset data, we can index the graph with an asset label:

graph[:ocgt]

GraphAssetData(:producer, true, true, 25.0, missing, 100.0, 0.0, 0.0, false, 0.0, 0.0, 0.0, 0.0, Dict{Symbol, Vector{Float64}}(), Dict{Tuple{Symbol, Int64}, Vector{Float64}}(), UnitRange{Int64}[], Dict{Int64, Vector{UnitRange{Int64}}}(2 => [1:1, 2:2, 3:3, 4:4, 5:5, 6:6, 7:7, 8:8, 9:9, 10:10  …  15:15, 16:16, 17:17, 18:18, 19:19, 20:20, 21:21, 22:22, 23:23, 24:24], 3 => [1:1, 2:2, 3:3, 4:4, 5:5, 6:6, 7:7, 8:8, 9:9, 10:10  …  15:15, 16:16, 17:17, 18:18, 19:19, 20:20, 21:21, 22:22, 23:23, 24:24], 1 => [1:1, 2:2, 3:3, 4:4, 5:5, 6:6, 7:7, 8:8, 9:9, 10:10  …  15:15, 16:16, 17:17, 18:18, 19:19, 20:20, 21:21, 22:22, 23:23, 24:24]), -1.0, Dict{Tuple{Int64, UnitRange{Int64}}, Float64}(), Dict{UnitRange{Int64}, Float64}())

This is a Julia struct, or composite type, named GraphAssetData. We can access its fields with .:

graph[:ocgt].type

:producer

Since labels returns a generator, we can iterate over its contents without collecting it into a vector.

for a in labels(graph)
    println("Asset $a has type $(graph[a].type)")
end

Asset ocgt has type producer
Asset ccgt has type producer
Asset wind has type producer
Asset solar has type producer
Asset ens has type producer
Asset demand has type consumer

To get all flows we can use edge_labels:

edge_labels(graph) |> collect

5-element Vector{Tuple{Symbol, Symbol}}:
 (:ocgt, :demand)
 (:ccgt, :demand)
 (:wind, :demand)
 (:solar, :demand)
 (:ens, :demand)

To access the flow data, we index with graph[u, v]:

graph[:ocgt, :demand]

GraphFlowData(:electricity, true, false, false, false, 0.07, 0.0, missing, 0.0, 0.0, 0.0, 0.0, Dict{Symbol, Vector{Float64}}(), Dict{Tuple{Symbol, Int64}, Vector{Float64}}(), UnitRange{Int64}[], Dict{Int64, Vector{UnitRange{Int64}}}(2 => [1:1, 2:2, 3:3, 4:4, 5:5, 6:6, 7:7, 8:8, 9:9, 10:10  …  15:15, 16:16, 17:17, 18:18, 19:19, 20:20, 21:21, 22:22, 23:23, 24:24], 3 => [1:1, 2:2, 3:3, 4:4, 5:5, 6:6, 7:7, 8:8, 9:9, 10:10  …  15:15, 16:16, 17:17, 18:18, 19:19, 20:20, 21:21, 22:22, 23:23, 24:24], 1 => [1:1, 2:2, 3:3, 4:4, 5:5, 6:6, 7:7, 8:8, 9:9, 10:10  …  15:15, 16:16, 17:17, 18:18, 19:19, 20:20, 21:21, 22:22, 23:23, 24:24]), Dict{Tuple{Int64, UnitRange{Int64}}, Float64}(), -1.0)

The type of the flow struct is GraphFlowData.

We can easily find all assets v for which a flow (a, v) exists for a given asset a (in this case, demand):

inneighbor_labels(graph, :demand) |> collect

5-element Vector{Symbol}:
 :ocgt
 :ccgt
 :wind
 :solar
 :ens

Similarly, all assets u for which a flow (u, a) exists for a given asset a (in this case, ocgt):

outneighbor_labels(graph, :ocgt) |> collect

1-element Vector{Symbol}:
 :demand

Manipulating the solution

First, see the description of the solution object.

Let's consider the larger dataset "Norse" in this section. And let's talk about two ways to access the solution.

The solution returned by solve_model

The solution, as shown before, can be obtained when calling solve_model or solve_model!.

using TulipaEnergyModel

# input_dir should be the path to Norse
energy_problem = create_energy_problem_from_csv_folder(input_dir)
create_model!(energy_problem)
solution = solve_model!(energy_problem)

To create a traditional array in the order given by the investable assets, one can run

using MetaGraphsNext

graph = energy_problem.graph
[solution.assets_investment[a] for a in labels(graph) if graph[a].investable]

12-element Vector{Float64}:
  1255.0
   500.0
    83.0
 14708.0
     1.7890356162559602e-5
     0.0
    80.0
    67.0
     0.0
   163.0
    33.0
     0.0

To create a traditional array in the order given by the investable flows, one can run

[solution.flows_investment[(u, v)] for (u, v) in edge_labels(graph) if graph[u, v].investable]

3-element Vector{Float64}:
 5.0
 8.510447821994182
 6.0

The solution.flow, solution.storage_level_intra_rp, and solution.storage_level_inter_rp values are linearized according to the dataframes in the dictionary energy_problem.dataframes with keys :flows, :lowest_storage_level_intra_rp, and :storage_level_inter_rp, respectively. You need to query the data from these dataframes and then use the column index to select the appropriate value.

To create a vector with all values of flow for a given (u, v) and rp, one can run

(u, v) = first(edge_labels(graph))
rp = 1
df = filter(
    row -> row.rp == rp && row.from == u && row.to == v,
    energy_problem.dataframes[:flows],
    view = true,
)
[solution.flow[row.index] for row in eachrow(df)]

168-element Vector{Float64}:
     0.0
     0.0
     0.0
     0.0
     0.0
     0.0
     0.0
     0.0
     0.0
     0.0
     ⋮
  5386.687839696075
  6583.403257693528
 10502.25248910553
  9977.343240986818
  7756.245816043869
  3186.1304959668196
     0.0
     0.0
     0.0

To create a vector with the all values of storage_level_intra_rp for a given a and rp, one can run

a = energy_problem.dataframes[:lowest_storage_level_intra_rp].asset[1]
rp = 1
df = filter(
    row -> row.asset == a && row.rp == rp,
    energy_problem.dataframes[:lowest_storage_level_intra_rp],
    view = true,
)
[solution.storage_level_intra_rp[row.index] for row in eachrow(df)]

168-element Vector{Float64}:
     0.0
     0.0
     0.0
     0.0
     0.0
     0.0
     0.0
    -0.0
  2089.9999999999995
  8075.0
     ⋮
 40005.65821031217
 33075.76004431899
 22020.7574242079
 11518.290854748093
  3353.8215747019153
     0.0
     0.0
     0.0
     0.0

To create a vector with the all values of storage_level_inter_rp for a given a, one can run

a = energy_problem.dataframes[:storage_level_inter_rp].asset[1]
df = filter(
    row -> row.asset == a,
    energy_problem.dataframes[:storage_level_inter_rp],
    view = true,
)
[solution.storage_level_inter_rp[row.index] for row in eachrow(df)]

52-element Vector{Float64}:
 25084.0
 25167.999999999996
 25251.999999999993
 25335.999999999993
 25419.99999999999
 25503.999999999985
 25587.999999999978
 25671.999999999975
 25755.99999999997
 25839.999999999967
     ⋮
 28978.947368421064
 29074.736842105278
 29158.736842105278
 29242.736842105278
 29326.736842105278
 29410.736842105278
 29494.736842105278
 29578.736842105278
 29674.736842105278

The solution inside the graph

In addition to the solution object, the solution is also stored by the individual assets and flows when solve_model! is called (i.e., when using an EnergyProblem object).

They can be accessed like any other value from GraphAssetData or GraphFlowData, which means that we recreate the values from the previous section in a new way:

[energy_problem.graph[a].investment for a in labels(graph) if graph[a].investable]

12-element Vector{Float64}:
  1255.0
   500.0
    83.0
 14708.0
     1.7890356162559602e-5
     0.0
    80.0
    67.0
     0.0
   163.0
    33.0
     0.0

[energy_problem.graph[u, v].investment for (u, v) in edge_labels(graph) if graph[u, v].investable]

3-element Vector{Float64}:
 5.0
 8.510447821994182
 6.0

(u, v) = first(edge_labels(graph))
rp = 1
df = filter(
    row -> row.rp == rp && row.from == u && row.to == v,
    energy_problem.dataframes[:flows],
    view = true,
)
[energy_problem.graph[u, v].flow[(rp, row.timesteps_block)] for row in eachrow(df)]

168-element Vector{Float64}:
     0.0
     0.0
     0.0
     0.0
     0.0
     0.0
     0.0
     0.0
     0.0
     0.0
     ⋮
  5386.687839696075
  6583.403257693528
 10502.25248910553
  9977.343240986818
  7756.245816043869
  3186.1304959668196
     0.0
     0.0
     0.0

To create a vector with all the values of storage_level_intra_rp for a given a and rp, one can run

a = energy_problem.dataframes[:lowest_storage_level_intra_rp].asset[1]
rp = 1
df = filter(
    row -> row.asset == a && row.rp == rp,
    energy_problem.dataframes[:lowest_storage_level_intra_rp],
    view = true,
)
[energy_problem.graph[a].storage_level_intra_rp[(rp, row.timesteps_block)] for row in eachrow(df)]

168-element Vector{Float64}:
     0.0
     0.0
     0.0
     0.0
     0.0
     0.0
     0.0
    -0.0
  2089.9999999999995
  8075.0
     ⋮
 40005.65821031217
 33075.76004431899
 22020.7574242079
 11518.290854748093
  3353.8215747019153
     0.0
     0.0
     0.0
     0.0

To create a vector with all the values of storage_level_inter_rp for a given a, one can run

a = energy_problem.dataframes[:storage_level_inter_rp].asset[1]
df = filter(
    row -> row.asset == a,
    energy_problem.dataframes[:storage_level_inter_rp],
    view = true,
)
[energy_problem.graph[a].storage_level_inter_rp[row.periods_block] for row in eachrow(df)]

52-element Vector{Float64}:
 25084.0
 25167.999999999996
 25251.999999999993
 25335.999999999993
 25419.99999999999
 25503.999999999985
 25587.999999999978
 25671.999999999975
 25755.99999999997
 25839.999999999967
     ⋮
 28978.947368421064
 29074.736842105278
 29158.736842105278
 29242.736842105278
 29326.736842105278
 29410.736842105278
 29494.736842105278
 29578.736842105278
 29674.736842105278

The solution inside the dataframes object

In addition to being stored in the solution object, and in the graph object, the solution for the flow, storage_level_intra_rp, and storage_level_inter_rp is also stored inside the corresponding DataFrame objects if solve_model! is called.

The code below will do the same as in the two previous examples:

(u, v) = first(edge_labels(graph))
rp = 1
df = filter(
    row -> row.rp == rp && row.from == u && row.to == v,
    energy_problem.dataframes[:flows],
    view = true,
)
df.solution

168-element view(::Vector{Float64}, 1:168) with eltype Float64:
     0.0
     0.0
     0.0
     0.0
     0.0
     0.0
     0.0
     0.0
     0.0
     0.0
     ⋮
  5386.687839696075
  6583.403257693528
 10502.25248910553
  9977.343240986818
  7756.245816043869
  3186.1304959668196
     0.0
     0.0
     0.0

a = energy_problem.dataframes[:storage_level_inter_rp].asset[1]
df = filter(
    row -> row.asset == a,
    energy_problem.dataframes[:storage_level_inter_rp],
    view = true,
)
df.solution

52-element view(::Vector{Float64}, 1:52) with eltype Float64:
 25084.0
 25167.999999999996
 25251.999999999993
 25335.999999999993
 25419.99999999999
 25503.999999999985
 25587.999999999978
 25671.999999999975
 25755.99999999997
 25839.999999999967
     ⋮
 28978.947368421064
 29074.736842105278
 29158.736842105278
 29242.736842105278
 29326.736842105278
 29410.736842105278
 29494.736842105278
 29578.736842105278
 29674.736842105278

a = energy_problem.dataframes[:lowest_storage_level_intra_rp].asset[1]
rp = 1
df = filter(
    row -> row.asset == a && row.rp == rp,
    energy_problem.dataframes[:lowest_storage_level_intra_rp],
    view = true,
)
df.solution

168-element view(::Vector{Float64}, 8:175) with eltype Float64:
     0.0
     0.0
     0.0
     0.0
     0.0
     0.0
     0.0
    -0.0
  2089.9999999999995
  8075.0
     ⋮
 40005.65821031217
 33075.76004431899
 22020.7574242079
 11518.290854748093
  3353.8215747019153
     0.0
     0.0
     0.0
     0.0

Values of constraints and expressions

By accessing the model directly, we can query the values of constraints and expressions. We need to know the name of the constraint and how it is indexed, and for that, you will need to check the model.

For instance, we can get all incoming flows in the lowest resolution for a given asset for a given representative period with the following:

using JuMP
a = energy_problem.dataframes[:lowest].asset[end]
rp = 1
df = filter(
    row -> row.asset == a && row.rp == rp,
    energy_problem.dataframes[:lowest],
    view = true,
)
[value(energy_problem.model[:incoming_flow_lowest_resolution][row.index]) for row in eachrow(df)]

168-element Vector{Float64}:
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 ⋮
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0

The values of constraints can also be obtained, however, they are frequently indexed in a subset, which means that their indexing is not straightforward. To know how they are indexed, it is necessary to look at the model code. For instance, to get the consumer balance, we first need to filter the :highest_in_out dataframes by consumers:

df_consumers = filter(
    row -> graph[row.asset].type == :consumer,
    energy_problem.dataframes[:highest_in_out],
    view = false,
);

We set view = false to create a copy of this DataFrame so we can make our indexes:

df_consumers.index = 1:size(df_consumers, 1) # overwrites existing index

1:960

Now we can filter this DataFrame. Note that the names in the stored dataframes are defined as Symbol.

a = :Asgard_E_demand
df = filter(
    row -> row.asset == a && row.rp == rp,
    df_consumers,
    view = true,
)
value.(energy_problem.model[:consumer_balance][df.index])

168-element Vector{Float64}:
 33951.9255800345
 32339.440666441926
 30618.08007858144
 29154.857719356212
 28466.28865613107
 27703.76592365483
 26184.454469129858
 26463.520980592206
 27178.204198305513
 29616.349838426668
     ⋮
 60971.61327059055
 61868.328688588
 65787.17792
 65262.26867188129
 63041.17124693834
 58471.05592686129
 54777.52475361511
 52128.801264625705
 46907.046645223265

Here value. (i.e., broadcasting) was used instead of the vector comprehension from previous examples just to show that it also works.

The value of the constraint is obtained by looking only at the part with variables. So a constraint like 2x + 3y - 1 <= 4 would return the value of 2x + 3y.

Writing the output to CSV

To save the solution to CSV files, you can use save_solution_to_file:

mkdir("outputs")
save_solution_to_file("outputs", energy_problem)

Plotting

In the previous sections, we have shown how to create vectors such as the one for flows. If you want simple plots, you can plot the vectors directly using any package you like.

If you would like more custom plots, check out TulipaPlots.jl, under development, which provides tailor-made plots for TulipaEnergyModel.jl.

Hydrothermal Dispatch example

Under development!