Fix problems with TwoLevelTree (#94)

* Minor simplification in checks
* Fixed problem with `scale_op_sp` and `TwoLevelTree`

Author: JulStraus

NEWS.md

@@ -1,5 +1,13 @@
 # Release notes
 
+## Version 0.10.4 (2026-06-15)
+
+### Bug fixes
+
+* Fixed a bug regarding the branch probability when using a `TwoLevelTree` structure:
+  * This bug resulted in strategic variables being scaled when calculated from operational variables through the function `scale_op_sp`.
+  * As a consequence, as an example, emission limits where wrongly applied.
+
 ## Version 0.10.3 (2026-06-10)
 
 ### Bug fixes

Project.toml

@@ -1,7 +1,7 @@
 name = "EnergyModelsBase"
 uuid = "5d7e687e-f956-46f3-9045-6f5a5fd49f50"
 authors = ["Lars Hellemo <Lars.Hellemo@sintef.no>, Julian Straus <Julian.Straus@sintef.no>"]
-version = "0.10.3"
+version = "0.10.4"
 
 [deps]
 JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
@@ -18,5 +18,5 @@ EMIExt = "EnergyModelsInvestments"
 EnergyModelsInvestments = "0.9"
 JuMP = "1"
 SparseVariables = "0.7.3"
-TimeStruct = "0.9"
+TimeStruct = "0.9.11"
 julia = "1.10"

docs/src/manual/optimization-variables.md

@@ -9,7 +9,7 @@ The latter is the recommended approach.
 !!! note
     The majority of the variables in `EnergyModelsBase` are rate variables.
     This imples that they are calculated for either an operational period duration of 1, when indexed over operational period ``t`` or a strategic period duration of 1, when indexed over strategic period ``t_\texttt{inv}``.
-    Typical units for rates are MW for energy streams, tonne/hour for mass streams, tonne/year for strategic emissions, and €/year for operational expenditures.
+    Typical units for rates are MW for energy streams, tonne/hour for mass streams, tonne/year for strategic emissions, and €/year for operating expenses.
     In this example, the duration of an operational period of 1 corresponds to an hour, while the duration of a strategic period of 1 corresponds to a year.
 
     Variables that are energy/mass based have that property highlighted in the documentation below.
@@ -27,9 +27,16 @@ The multiplication then leads to an energy/mass quantity in stead of an energy/m
 The coupling of strategic and operational periods can be achieved through the function `scale_op_sp(t, t_inv)`.
 This functions allows for considering the scaling of the operational periods within a strategic period.
 
+!!! note "`TwoLevelTree` and variables"
+    All variables that are indexed over strategic periods do not take into consideration the branch probability.
+    This is, *e.g.*, the case for the strategic emission variables or the variable operating expenses that apply per scenario/branch.
+    The reason for this approach is to simplify the comparison of individual values between different branches without the need to consider the probability.
+
+    The branch probability is however taken into account when calculating the objective function.
+
 ## [Operational cost variables](@id man-opt_var-opex)
 
-Operational cost variables are included to account for operational expenditures (OPEX) of the model.
+Operational cost variables are included to account for operating expenses (OPEX) of the model.
 These costs are pure dependent on either the use or the installed capacity of a node ``n``.
 All nodes ``n`` (except [`Availability`](@ref)-nodes) have the following variables representing the operational costs of the nodes:
 

src/checks.jl

@@ -710,12 +710,12 @@ function check_representative_profile(time_profile::TimeProfile, message::String
     # Iterate through the strategic profiles, if existing
     if isa(time_profile, StrategicProfile)
         for l1_profile ∈ time_profile.vals
-            sub_msg = "in strategic profiles " * message
-            bool_rp = check_repr_sub_profile(l1_profile, sub_msg, bool_rp)
+            sub_msg_1 = "in strategic profiles " * message
+            bool_rp = check_repr_sub_profile(l1_profile, sub_msg_1, bool_rp)
             if isa(l1_profile, RepresentativeProfile)
                 for l2_profile ∈ l1_profile.vals
-                    sub_msg = "in representative profiles in strategic profiles " * message
-                    bool_rp = check_repr_sub_profile(l2_profile, sub_msg, bool_rp)
+                    sub_msg_2 = "in representative profiles " * sub_msg_1
+                    bool_rp = check_repr_sub_profile(l2_profile, sub_msg_2, bool_rp)
                 end
             end
         end
@@ -748,8 +748,8 @@ Function for checking that an individual `TimeProfile` does not include the wron
 scenario indexing.
 
 ## Checks
-- `TimeProfile`s accessed in `RepresentativePeriod`s cannot include `OperationalProfile`
-  or `ScenarioProfile` as this is not allowed through indexing on the `TimeProfile`.
+- `TimeProfile`s accessed in `OperationalScenario`s cannot include `OperationalProfile` as
+   this is not allowed through indexing on the `TimeProfile`.
 """
 function check_scenario_profile(time_profile::TimeProfile, message::String)
     # Check on the highest level
@@ -758,23 +758,23 @@ function check_scenario_profile(time_profile::TimeProfile, message::String)
     # Iterate through the strategic profiles, if existing
     if isa(time_profile, StrategicProfile)
         for l1_profile ∈ time_profile.vals
-            sub_msg = "in strategic profiles " * message
-            bool_scp = check_osc_sub_profile(l1_profile, sub_msg, bool_scp)
+            sub_msg_1 = "in strategic profiles " * message
+            bool_scp = check_osc_sub_profile(l1_profile, sub_msg_1, bool_scp)
             if isa(l1_profile, RepresentativeProfile)
+                sub_msg_2 = "in representative profiles " * sub_msg_1
                 for l2_profile ∈ l1_profile.vals
-                    sub_msg = "in representative profiles in strategic profiles " * message
-                    bool_scp = check_osc_sub_profile(l2_profile, sub_msg, bool_scp)
+                    bool_scp = check_osc_sub_profile(l2_profile, sub_msg_2, bool_scp)
                     if isa(l2_profile, ScenarioProfile)
+                        sub_msg_3 = "in scenario profiles in " * sub_msg_2
                         for l3_profile ∈ l2_profile.vals
-                            sub_msg = "in scenario profiles in representative profiles in strategic profiles " * message
-                            bool_scp = check_osc_sub_profile(l3_profile, sub_msg, bool_scp)
+                            bool_scp = check_osc_sub_profile(l3_profile, sub_msg_3, bool_scp)
                         end
                     end
                 end
             elseif isa(l1_profile, ScenarioProfile)
                 for l2_profile ∈ l1_profile.vals
-                    sub_msg = "in scenario profiles in strategic profiles " * message
-                    bool_scp = check_osc_sub_profile(l2_profile, sub_msg, bool_scp)
+                    sub_msg_2 = "in scenario profiles " * sub_msg_1
+                    bool_scp = check_osc_sub_profile(l2_profile, sub_msg_2, bool_scp)
                 end
             end
         end
@@ -783,12 +783,12 @@ function check_scenario_profile(time_profile::TimeProfile, message::String)
     # Iterate through the representative profiles, if existing
     if isa(time_profile, RepresentativeProfile)
         for l1_profile ∈ time_profile.vals
-            sub_msg = "in representative profiles " * message
-            bool_scp = check_osc_sub_profile(l1_profile, sub_msg, bool_scp)
+            sub_msg_1 = "in representative profiles " * message
+            bool_scp = check_osc_sub_profile(l1_profile, sub_msg_1, bool_scp)
             if isa(l1_profile, ScenarioProfile)
                 for l2_profile ∈ l1_profile.vals
-                    sub_msg = "in scenario profiles in representative profiles " * message
-                    bool_scp = check_osc_sub_profile(l2_profile, sub_msg, bool_scp)
+                    sub_msg_2 = "in scenario profiles " * sub_msg_1
+                    bool_scp = check_osc_sub_profile(l2_profile, sub_msg_2, bool_scp)
                 end
             end
         end

src/model.jl

@@ -762,7 +762,7 @@ function objective(m, 𝒳ᵛᵉᶜ, 𝒫, 𝒯, modeltype::EnergyModel)
     # Calculation of the objective function.
     @objective(m, Max,
         -sum(
-            sum(𝒳[t_inv] for 𝒳 ∈ opex) * duration_strat(t_inv)
+            sum(𝒳[t_inv] for 𝒳 ∈ opex) * duration_strat(t_inv) * probability_branch(t_inv)
         for t_inv ∈ 𝒯ᴵⁿᵛ)
     )
 end

src/utils.jl

@@ -319,11 +319,23 @@ end
 
 Provides a simplified function for returning the multiplication
 
-``duration(t) * multiple\\_strat(t\\_inv, t) * probability(t)``
+``duration(t) * multiple\\_strat(t\\_inv, t) * probability(t) / probability_branch(t)``
 
 when operational periods are coupled with strategic periods. It is used to scale the value
 provided for operational periods to a duration of 1 of a strategic period.
 
+!!! note "`TwoLevelTree` application"
+    The function does not consider the probability of a branch when using a
+    [`TwoLevelTree`](@extref TimeStruct.TwoLevelTree) time structure. The reason is that we
+    do not consider any scaling or discounting for the individual strategic variables to
+    avoid scaling the corresponding bounds. In addition, it allows a simple comparison
+    of branches with different probabilities.
+
+    The scaling is however included in the function [`objective`](@ref) by utilizing the
+    function [`probability_branch`](@extref TimeStruct.probability_branch)
+    or [`objective_weight`](@extref TimeStruct.objective_weight) when considering
+    investments.
+
 # Example
 
 ```julia
@@ -340,7 +352,7 @@ scale_op_sp(t_inv, t)
 ```
 """
 scale_op_sp(t_inv::TS.AbstractStrategicPeriod, t::TS.TimePeriod) =
-    duration(t) * multiple_strat(t_inv, t) * probability(t)
+    duration(t) * multiple_strat(t_inv, t) * probability(t) / probability_branch(t)
 
 function multiple(t_inv, t)
     @warn(

test/test_general.jl

@@ -1,22 +1,22 @@
-function generate_data()
+function generate_data(; 𝒯 = TwoLevel(4, 2, SimpleTimes(4, 2), op_per_strat = 8.0))
 
     # Define the different resources
     NG = ResourceEmit("NG", 0.2)
     Coal = ResourceCarrier("Coal", 0.35)
     Power = ResourceCarrier("Power", 0.0)
     CO2 = ResourceEmit("CO2", 1.0)
-    products = [NG, Coal, Power, CO2]
+    𝒫 = [NG, Coal, Power, CO2]
 
     # Creation of the emission data for the individual nodes.
     capture_data = CaptureEnergyEmissions(0.9)
     emission_data = EmissionsEnergy()
 
     # Create the individual test nodes, corresponding to a system with an electricity demand/sink,
     # coal and nautral gas sources, coal and natural gas (with CCS) power plants and CO2 storage.
-    nodes = [
-        GenAvailability(1, products),
-        RefSource(2, FixedProfile(1e12), FixedProfile(30), FixedProfile(0), Dict(NG => 1)),
-        RefSource(3, FixedProfile(1e12), FixedProfile(9), FixedProfile(0), Dict(Coal => 1)),
+    𝒩 = [
+        GenAvailability(1, 𝒫),
+        RefSource(2, FixedProfile(100), FixedProfile(30), FixedProfile(0), Dict(NG => 1)),
+        RefSource(3, FixedProfile(100), FixedProfile(9), FixedProfile(0), Dict(Coal => 1)),
         RefNetworkNode(
             4,
             FixedProfile(25),
@@ -52,34 +52,33 @@ function generate_data()
     ]
 
     # Connect all nodes with the availability node for the overall energy/mass balance
-    links = [
-        Direct(14, nodes[1], nodes[4], Linear())
-        Direct(15, nodes[1], nodes[5], Linear())
-        Direct(16, nodes[1], nodes[6], Linear())
-        Direct(17, nodes[1], nodes[7], Linear())
-        Direct(21, nodes[2], nodes[1], Linear())
-        Direct(31, nodes[3], nodes[1], Linear())
-        Direct(41, nodes[4], nodes[1], Linear())
-        Direct(51, nodes[5], nodes[1], Linear())
-        Direct(61, nodes[6], nodes[1], Linear())
+    ℒ = [
+        Direct(14, 𝒩[1], 𝒩[4], Linear())
+        Direct(15, 𝒩[1], 𝒩[5], Linear())
+        Direct(16, 𝒩[1], 𝒩[6], Linear())
+        Direct(17, 𝒩[1], 𝒩[7], Linear())
+        Direct(21, 𝒩[2], 𝒩[1], Linear())
+        Direct(31, 𝒩[3], 𝒩[1], Linear())
+        Direct(41, 𝒩[4], 𝒩[1], Linear())
+        Direct(51, 𝒩[5], 𝒩[1], Linear())
+        Direct(61, 𝒩[6], 𝒩[1], Linear())
     ]
 
-    # Creation of the time structure and global data
-    T = TwoLevel(4, 2, SimpleTimes(4, 2), op_per_strat = 8)
-    model = OperationalModel(
+    # Creation of the modeltype
+    modeltype = OperationalModel(
         Dict(CO2 => StrategicProfile([160, 140, 120, 100]), NG => FixedProfile(1e6)),
         Dict(CO2 => FixedProfile(10)),
         CO2,
     )
 
     # Input data structure
-    case = Case(T, products, [nodes, links], [[get_nodes, get_links]])
-    return case, model
+    case = Case(𝒯, 𝒫, [𝒩, ℒ], [[get_nodes, get_links]])
+    return case, modeltype
 end
 
 @testset "General tests" begin
-    case, model = generate_data()
-    m = run_model(case, model, HiGHS.Optimizer)
+    case, modeltype = generate_data()
+    m = run_model(case, modeltype, HiGHS.Optimizer)
 
     # Retrieve data from the case structure
     𝒫 = get_products(case)
@@ -100,6 +99,8 @@ end
 
     ℒ = get_links(case)
 
+    objective_TwoLevel = objective_value(m)
+
     # Check for the objective value
     # (*2 compared to 0.6.0 due to change in strategic period duration)
     # (-10400 = 2*10*(160+140+120+100) compared to 0.8.3 due to inclusion of co2 emissions)
@@ -113,10 +114,10 @@ end
     # - constraints_emissions(m, 𝒩, 𝒯, 𝒫, modeltype::EnergyModel)
     @test all(
         value.(m[:emissions_strategic])[t_inv, CO2] <=
-        EMB.emission_limit(model, CO2, t_inv) for t_inv ∈ 𝒯ᴵⁿᵛ
+        EMB.emission_limit(modeltype, CO2, t_inv) for t_inv ∈ 𝒯ᴵⁿᵛ
     )
     @test all(
-        value.(m[:emissions_strategic])[t_inv, NG] <= EMB.emission_limit(model, NG, t_inv)
+        value.(m[:emissions_strategic])[t_inv, NG] <= EMB.emission_limit(modeltype, NG, t_inv)
         for t_inv ∈ 𝒯ᴵⁿᵛ
     )
 
@@ -140,7 +141,7 @@ end
                 value.(m[:opex_var][n, t_inv]) + value.(m[:opex_fixed][n, t_inv])
             for n ∈ 𝒩ᵒᵖᵉˣ) +
             sum(
-                value.(m[:emissions_total][t, CO2]) * emission_price(model, CO2, t) *
+                value.(m[:emissions_total][t, CO2]) * emission_price(modeltype, CO2, t) *
                 scale_op_sp(t_inv, t) for t ∈ t_inv)
         ) * duration_strat(t_inv) for t_inv ∈ 𝒯ᴵⁿᵛ
     ) ≈ objective_value(m) atol = TEST_ATOL
@@ -184,4 +185,10 @@ end
             p ∈ EMB.link_res(l), atol ∈ TEST_ATOL
         ) for l ∈ ℒ, atol ∈ TEST_ATOL
     )
+
+    # Test that the results are exactly the same for an equivalent `TwoLevelTree`
+    𝒯 = TwoLevelTree(2, [2, 2, 2], SimpleTimes(4, 2), op_per_strat = 8.0)
+    case, modeltype = generate_data(; 𝒯)
+    m = run_model(case, modeltype, HiGHS.Optimizer)
+    @test objective_TwoLevel ≈ objective_value(m)
 end