[WIP] Fix the duration function and the storage constraints

src/model.jl

@@ -19,7 +19,7 @@ function create_model(graph, params, sets; verbose = false, write_lp_file = fals
     Fi  = [f for f ∈ F if params.flows_investable[f]]
     Ft  = [f for f ∈ F if params.flows_is_transport[f]]
     # K_rp = sets.time_steps
-    K_A = sets.time_intervals_per_asset
+    # K_A = sets.time_intervals_per_asset
     K_F = sets.time_intervals_per_flow
     P = sets.constraints_time_periods
     RP = sets.rep_periods
@@ -32,7 +32,7 @@ function create_model(graph, params, sets; verbose = false, write_lp_file = fals
     @variable(model, flow[f ∈ F, rp ∈ RP, K_F[(f, rp)]])         #flow from asset a to asset aa [MW]
     @variable(model, 0 ≤ assets_investment[Ai], Int)  #number of installed asset units [N]
     @variable(model, 0 ≤ flows_investment[Fi], Int)
-    @variable(model, 0 ≤ storage_level[a ∈ As, rp ∈ RP, K_A[(a, rp)]])
+    @variable(model, 0 ≤ storage_level[a ∈ As, rp ∈ RP, P[(a, rp)]])
 
     # TODO: Fix storage_level[As, RP, 0] = 0
 
@@ -71,11 +71,11 @@ function create_model(graph, params, sets; verbose = false, write_lp_file = fals
     )
 
     # Constraints
-    # Computes the percentage of the `interval` that is within the `period`, and multiply by the
+    # Computes the duration of the `interval` that is within the `period`, and multiply by the
     # scale of the representative period `rp`.
     # It is equivalent to finding the indexes of these values in the matrix.
-    function ratio(T, I, rp)
-        return length(T ∩ I) / length(I) * params.time_scale[rp]
+    function duration(T, I, rp)
+        return length(T ∩ I) * params.time_scale[rp]
     end
 
     # Balance equations
@@ -84,10 +84,10 @@ function create_model(graph, params, sets; verbose = false, write_lp_file = fals
         model,
         c_consumer_balance[a ∈ Ac, rp ∈ RP, T ∈ P[(a, rp)]],
         sum(
-            ratio(T, I, rp) * flow[f, rp, I] for
+            duration(T, I, rp) * flow[f, rp, I] for
             f ∈ F, I ∈ sets.time_intervals_per_flow[(f, rp)] if f[2] == a
         ) - sum(
-            ratio(T, I, rp) * flow[f, rp, I] for
+            duration(T, I, rp) * flow[f, rp, I] for
             f ∈ F, I ∈ sets.time_intervals_per_flow[(f, rp)] if f[1] == a
         ) ==
         sum(get(params.assets_profile, (a, rp, k), 1.0) for k ∈ T) * params.peak_demand[a]
@@ -97,19 +97,13 @@ function create_model(graph, params, sets; verbose = false, write_lp_file = fals
     # TODO: Add p^{inflow}
     @constraint(
         model,
-        c_storage_balance[a ∈ As, rp ∈ RP, T ∈ P[(a, rp)]],
-        sum(
-            ratio([T[end]], I, rp) * storage_level[a, rp, I] for
-            I ∈ sets.time_intervals_per_asset[(a, rp)]
-        ) ==
-        sum(
-            ratio([T[1] - 1], I, rp) * storage_level[a, rp, I] for
-            I ∈ sets.time_intervals_per_asset[(a, rp)]
-        ) + sum(
-            ratio(T, I, rp) * flow[f, rp, I] * params.flows_efficiency[f] for
+        c_storage_balance[a ∈ As, rp ∈ RP, (k, T) ∈ enumerate(P[(a, rp)])],
+        storage_level[a, rp, T] ==
+        (k > 1 ? storage_level[a, rp, P[(a, rp)][k-1]] : 0.0) + sum(
+            duration(T, I, rp) * flow[f, rp, I] * params.flows_efficiency[f] for
             f ∈ F, I ∈ sets.time_intervals_per_flow[(f, rp)] if f[2] == a
         ) - sum(
-            ratio(T, I, rp) * flow[f, rp, I] / params.flows_efficiency[f] for
+            duration(T, I, rp) * flow[f, rp, I] / params.flows_efficiency[f] for
             f ∈ F, I ∈ sets.time_intervals_per_flow[(f, rp)] if f[1] == a
         )
     )
@@ -119,10 +113,10 @@ function create_model(graph, params, sets; verbose = false, write_lp_file = fals
         model,
         c_hub_balance[a ∈ Ah, rp ∈ RP, T ∈ P[(a, rp)]],
         sum(
-            ratio(T, I, rp) * flow[f, rp, I] for
+            duration(T, I, rp) * flow[f, rp, I] for
             f ∈ F, I ∈ sets.time_intervals_per_flow[(f, rp)] if f[2] == a
         ) == sum(
-            ratio(T, I, rp) * flow[f, rp, I] for
+            duration(T, I, rp) * flow[f, rp, I] for
             f ∈ F, I ∈ sets.time_intervals_per_flow[(f, rp)] if f[1] == a
         )
     )
@@ -132,10 +126,10 @@ function create_model(graph, params, sets; verbose = false, write_lp_file = fals
         model,
         c_conversion_balance[a ∈ Acv, rp ∈ RP, T ∈ P[(a, rp)]],
         sum(
-            ratio(T, I, rp) * flow[f, rp, I] * params.flows_efficiency[f] for
+            duration(T, I, rp) * flow[f, rp, I] * params.flows_efficiency[f] for
             f ∈ F, I ∈ sets.time_intervals_per_flow[(f, rp)] if f[2] == a
         ) == sum(
-            ratio(T, I, rp) * flow[f, rp, I] / params.flows_efficiency[f] for
+            duration(T, I, rp) * flow[f, rp, I] / params.flows_efficiency[f] for
             f ∈ F, I ∈ sets.time_intervals_per_flow[(f, rp)] if f[1] == a
         )
     )
@@ -146,7 +140,7 @@ function create_model(graph, params, sets; verbose = false, write_lp_file = fals
         model,
         c_overall_output_flows[a ∈ Acv∪As∪Ap, rp ∈ RP, T ∈ P[(a, rp)]],
         sum(
-            ratio(T, I, rp) * flow[f, rp, I] for
+            duration(T, I, rp) * flow[f, rp, I] for
             f ∈ F, I ∈ sets.time_intervals_per_flow[(f, rp)] if f[1] == a
         ) ≤
         sum(get(params.assets_profile, (a, rp, k), 1.0) for k ∈ T) * (
@@ -160,7 +154,7 @@ function create_model(graph, params, sets; verbose = false, write_lp_file = fals
         model,
         c_overall_input_flows[a ∈ As, rp ∈ RP, T ∈ P[(a, rp)]],
         sum(
-            ratio(T, I, rp) * flow[f, rp, I] for
+            duration(T, I, rp) * flow[f, rp, I] for
             f ∈ F, I ∈ sets.time_intervals_per_flow[(f, rp)] if f[2] == a
         ) ≤
         sum(get(params.assets_profile, (a, rp, k), 1.0) for k ∈ T) * (
@@ -180,7 +174,8 @@ function create_model(graph, params, sets; verbose = false, write_lp_file = fals
             !(a ∈ Ah ∪ Ac) && f[1] == a && f ∉ Ft,
         ],
         sum(
-            ratio(T, I, rp) * flow[f, rp, I] for I ∈ sets.time_intervals_per_flow[(f, rp)]
+            duration(T, I, rp) * flow[f, rp, I] for
+            I ∈ sets.time_intervals_per_flow[(f, rp)]
         ) ≤
         sum(get(params.assets_profile, (a, rp, k), 1.0) for k ∈ T) * (
             params.assets_init_capacity[a] +
@@ -235,10 +230,8 @@ function create_model(graph, params, sets; verbose = false, write_lp_file = fals
     @constraint(
         model,
         upper_bound_for_storage_level[a ∈ As, rp ∈ RP, T ∈ P[(a, rp)]],
-        sum(
-            ratio(T, I, rp) * storage_level[a, rp, I] for
-            I ∈ sets.time_intervals_per_asset[(a, rp)]
-        ) ≤ params.initial_storage_capacity[a] + (a ∈ Ai ? energy_limit[a] : 0.0)
+        storage_level[a, rp, T] ≤
+        params.initial_storage_capacity[a] + (a ∈ Ai ? energy_limit[a] : 0.0)
     )
 
     if write_lp_file