Compute mean realized flows

core/src/allocation_optim.jl

@@ -203,7 +203,7 @@ function set_initial_capacities_source!(
 )::Nothing
     (; problem) = allocation_model
     (; graph, allocation) = p
-    (; mean_flows) = allocation
+    (; mean_input_flows) = allocation
     (; subnetwork_id) = allocation_model
     source_constraints = problem[:source]
     main_network_source_edges = get_main_network_connections(p, subnetwork_id)
@@ -214,7 +214,7 @@ function set_initial_capacities_source!(
             # If it is a source edge for this allocation problem
             if edge ∉ main_network_source_edges
                 # Reset the source to the averaged flow over the last allocation period
-                source_capacity = mean_flows[edge][]
+                source_capacity = mean_input_flows[edge][]
                 JuMP.set_normalized_rhs(
                     source_constraints[edge],
                     # It is assumed that the allocation procedure does not have to be differentiated.
@@ -309,11 +309,11 @@ function get_basin_data(
     (; graph, basin, level_demand, allocation) = p
     (; vertical_flux) = basin
     (; Δt_allocation) = allocation_model
-    (; mean_flows) = allocation
+    (; mean_input_flows) = allocation
     @assert node_id.type == NodeType.Basin
     vertical_flux = get_tmp(vertical_flux, 0)
     _, basin_idx = id_index(basin.node_id, node_id)
-    influx = mean_flows[(node_id, node_id)][]
+    influx = mean_input_flows[(node_id, node_id)][]
     _, basin_idx = id_index(basin.node_id, node_id)
     storage_basin = u.storage[basin_idx]
     control_inneighbors = inneighbor_labels_type(graph, node_id, EdgeType.control)
@@ -709,7 +709,7 @@ function save_demands_and_allocations!(
     priority_idx::Int,
 )::Nothing
     (; graph, allocation, user_demand, flow_demand, basin) = p
-    (; record_demand, priorities) = allocation
+    (; record_demand, priorities, mean_realized_flows) = allocation
     (; subnetwork_id, problem) = allocation_model
     node_ids = graph[].node_ids[subnetwork_id]
     constraints_outflow = problem[:basin_outflow]
@@ -725,8 +725,7 @@ function save_demands_and_allocations!(
             user_demand_idx = findsorted(user_demand.node_id, node_id)
             demand = user_demand.demand[user_demand_idx, priority_idx]
             allocated = user_demand.allocated[user_demand_idx, priority_idx]
-            #NOTE: instantaneous
-            realized = get_flow(graph, inflow_id(graph, node_id), node_id, 0)
+            realized = mean_realized_flows[(inflow_id(graph, node_id), node_id)]
 
         elseif has_external_demand(graph, node_id, :level_demand)[1]
             basin_priority_idx = get_external_priority_idx(p, node_id)
@@ -745,9 +744,7 @@ function save_demands_and_allocations!(
                 end
                 allocated =
                     JuMP.value(F_basin_in[node_id]) - JuMP.value(F_basin_out[node_id])
-                # TODO: realized for a basin is not so clear, maybe it should be Δstorage/Δt
-                # over the last allocation interval?
-                realized = 0.0
+                realized = mean_realized_flows[(node_id, node_id)]
             end
 
         else
@@ -763,8 +760,7 @@ function save_demands_and_allocations!(
                         flow_demand_node_id,
                     )] : 0.0
                 allocated = JuMP.value(F[(inflow_id(graph, node_id), node_id)])
-                #NOTE: Still instantaneous
-                realized = get_flow(graph, inflow_id(graph, node_id), node_id, 0)
+                realized = mean_realized_flows[(inflow_id(graph, node_id), node_id)]
             end
         end
 
@@ -777,9 +773,6 @@ function save_demands_and_allocations!(
             push!(record_demand.priority, priorities[priority_idx])
             push!(record_demand.demand, demand)
             push!(record_demand.allocated, allocated)
-
-            # TODO: This is now the last abstraction before the allocation update,
-            # should be the average abstraction since the last allocation solve
             push!(record_demand.realized, realized)
         end
     end

core/src/callback.jl

@@ -114,21 +114,32 @@ function integrate_flows!(u, t, integrator)::Nothing
     end
 
     # Allocation source flows
-    for (edge, value) in allocation.mean_flows
+    for (edge, value) in allocation.mean_input_flows
         if edge[1] == edge[2]
             # Vertical fluxes
             _, basin_idx = id_index(basin.node_id, edge[1])
-            value[] +=
+            allocation.mean_input_flows[edge] =
+                value +
                 0.5 *
                 (get_influx(basin, basin_idx) + get_influx(basin, basin_idx; prev = true)) *
                 dt
         else
             # Horizontal flows
-            value[] +=
+            allocation.mean_input_flows[edge] =
+                value +
                 0.5 * (get_flow(graph, edge..., 0) + get_flow_prev(graph, edge..., 0)) * dt
         end
     end
 
+    # Realized demand flows
+    for (edge, value) in allocation.mean_realized_flows
+        if edge[1] !== edge[2]
+            value +=
+                0.5 * (get_flow(graph, edge..., 0) + get_flow_prev(graph, edge..., 0)) * dt
+            allocation.mean_realized_flows[edge] = value
+        end
+    end
+
     copyto!(flow_prev, flow)
     copyto!(vertical_flux_prev, vertical_flux)
     return nothing
@@ -460,8 +471,8 @@ end
 "Solve the allocation problem for all demands and assign allocated abstractions."
 function update_allocation!(integrator)::Nothing
     (; p, t, u) = integrator
-    (; allocation) = p
-    (; allocation_models, mean_flows) = allocation
+    (; allocation, basin) = p
+    (; allocation_models, mean_input_flows, mean_realized_flows) = allocation
 
     # Don't run the allocation algorithm if allocation is not active
     # (Specifically for running Ribasim via the BMI)
@@ -471,10 +482,21 @@ function update_allocation!(integrator)::Nothing
 
     (; Δt_allocation) = allocation_models[1]
 
-    # Divide by the allocation Δt to obtain the mean flows
+    # Divide by the allocation Δt to obtain the mean input flows
     # from the integrated flows
-    for value in values(mean_flows)
-        value[] /= Δt_allocation
+    for key in keys(mean_input_flows)
+        mean_input_flows[key] /= Δt_allocation
+    end
+
+    # Divide by the allocation Δt to obtain the mean realized flows
+    # from the integrated flows
+    for (edge, value) in mean_realized_flows
+        if edge[1] == edge[2]
+            # Compute the mean realized demand for basins as Δstorage/Δt_allocation
+            _, basin_idx = id_index(basin.node_id, edge[1])
+            mean_realized_flows[edge] = value + u[basin_idx]
+        end
+        mean_realized_flows[edge] /= Δt_allocation
     end
 
     # If a main network is present, collect demands of subnetworks
@@ -491,9 +513,19 @@ function update_allocation!(integrator)::Nothing
         allocate_demands!(p, allocation_model, t, u)
     end
 
-    # Reset the mean source flows
-    for value in values(mean_flows)
-        value[] = 0.0
+    # Reset the mean flows
+    for mean_flows in (mean_input_flows, mean_realized_flows)
+        for edge in keys(mean_flows)
+            mean_flows[edge] = 0.0
+        end
+    end
+
+    # Set basin storages for mean storage change computation
+    for (edge, value) in mean_realized_flows
+        if edge[1] == edge[2]
+            _, basin_idx = id_index(basin.node_id, edge[1])
+            mean_realized_flows[edge] = value - u[basin_idx]
+        end
     end
 end
 

core/src/parameter.jl

@@ -63,7 +63,8 @@ main_network_connections: (from_id, to_id) from the main network to the subnetwo
 priorities: All used priority values.
 subnetwork_demands: The demand of an edge from the main network to a subnetwork
 subnetwork_allocateds: The allocated flow of an edge from the main network to a subnetwork
-mean_flows: Flows averaged over Δt_allocation over edges that are allocation sources
+mean_input_flows: Flows averaged over Δt_allocation over edges that are allocation sources
+mean_realized_flows: Flows averaged over Δt_allocation over edges that realize a demand
 record_demand: A record of demands and allocated flows for nodes that have these
 record_flow: A record of all flows computed by allocation optimization, eventually saved to
     output file
@@ -75,7 +76,8 @@ struct Allocation
     priorities::Vector{Int32}
     subnetwork_demands::Dict{Tuple{NodeID, NodeID}, Vector{Float64}}
     subnetwork_allocateds::Dict{Tuple{NodeID, NodeID}, Vector{Float64}}
-    mean_flows::Dict{Tuple{NodeID, NodeID}, Base.RefValue{Float64}}
+    mean_input_flows::Dict{Tuple{NodeID, NodeID}, Float64}
+    mean_realized_flows::Dict{Tuple{NodeID, NodeID}, Float64}
     record_demand::@NamedTuple{
         time::Vector{Float64},
         subnetwork_id::Vector{Int32},

core/src/read.jl

@@ -1036,20 +1036,40 @@ function Allocation(db::DB, config::Config, graph::MetaGraph)::Allocation
         optimization_type = String[],
     )
 
-    mean_flows = Dict{Tuple{NodeID, NodeID}, Base.RefValue{Float64}}()
+    mean_input_flows = Dict{Tuple{NodeID, NodeID}, Float64}()
 
     # Find edges which serve as sources in allocation
     for edge_metadata in values(graph.edge_data)
         (; subnetwork_id_source, edge) = edge_metadata
         if subnetwork_id_source != 0
-            mean_flows[edge] = Ref(0.0)
+            mean_input_flows[edge] = 0.0
         end
     end
 
     # Find basins with a level demand
     for node_id in values(graph.vertex_labels)
         if has_external_demand(graph, node_id, :level_demand)[1]
-            mean_flows[(node_id, node_id)] = Ref(0.0)
+            mean_input_flows[(node_id, node_id)] = 0.0
+        end
+    end
+
+    mean_realized_flows = Dict{Tuple{NodeID, NodeID}, Float64}()
+
+    # Find edges that realize a demand
+    for edge_metadata in values(graph.edge_data)
+        (; type, edge) = edge_metadata
+
+        src_id, dst_id = edge
+        user_demand_inflow = (type == EdgeType.flow) && (dst_id.type == NodeType.UserDemand)
+        level_demand_inflow =
+            (type == EdgeType.control) && (src_id.type == NodeType.LevelDemand)
+        flow_demand_inflow =
+            (type == EdgeType.flow) && has_external_demand(graph, dst_id, :flow_demand)[1]
+
+        if user_demand_inflow || flow_demand_inflow
+            mean_realized_flows[edge] = 0.0
+        elseif level_demand_inflow
+            mean_realized_flows[(dst_id, dst_id)] = 0.0
         end
     end
 
@@ -1060,7 +1080,8 @@ function Allocation(db::DB, config::Config, graph::MetaGraph)::Allocation
         get_all_priorities(db, config),
         Dict{Tuple{NodeID, NodeID}, Vector{Float64}}(),
         Dict{Tuple{NodeID, NodeID}, Vector{Float64}}(),
-        mean_flows,
+        mean_input_flows,
+        mean_realized_flows,
         record_demand,
         record_flow,
     )

core/src/util.jl

@@ -758,7 +758,7 @@ as input. Therefore we set the instantaneous flows as the mean flows as allocati
 function set_initial_allocation_mean_flows!(integrator)::Nothing
     (; u, p, t) = integrator
     (; allocation, graph, basin) = p
-    (; mean_flows, allocation_models) = allocation
+    (; mean_input_flows, mean_realized_flows, allocation_models) = allocation
     (; Δt_allocation) = allocation_models[1]
     (; vertical_flux) = basin
     vertical_flux = get_tmp(vertical_flux, 0)
@@ -768,15 +768,28 @@ function set_initial_allocation_mean_flows!(integrator)::Nothing
     # one is just to make sure.
     water_balance!(get_du(integrator), u, p, t)
 
-    for (edge, value) in mean_flows
+    for edge in keys(mean_input_flows)
         if edge[1] == edge[2]
             q = get_influx(basin, edge[1])
         else
             q = get_flow(graph, edge..., 0)
         end
         # Multiply by Δt_allocation as averaging divides by this factor
         # in update_allocation!
-        value[] = q * Δt_allocation
+        mean_input_flows[edge] = q * Δt_allocation
     end
+
+    # Mean realized demands for basins are calculated as Δstorage/Δt
+    # This sets the realized demands as -storage_old
+    for edge in keys(mean_realized_flows)
+        if edge[1] == edge[2]
+            _, basin_idx = id_index(basin.node_id, edge[1])
+            mean_realized_flows[edge] = -u[basin_idx]
+        else
+            q = get_flow(graph, edge..., 0)
+            mean_realized_flows[edge] = q * Δt_allocation
+        end
+    end
+
     return nothing
 end