Improved testing for subnetwork communication

core/src/allocation_init.jl

@@ -945,7 +945,7 @@ function has_demand_priority_subnetwork(
     p_independent::ParametersIndependent,
     node_ids_in_subnetwork::NodeIDsInSubnetwork,
 )::Vector{Bool}
-    (; allocation, graph, user_demand, flow_demand, level_demand) = p_independent
+    (; allocation, user_demand, flow_demand, level_demand) = p_independent
     (; demand_priorities_all) = allocation
     (;
         user_demand_ids_subnetwork,

core/src/allocation_optim.jl

@@ -345,7 +345,10 @@ function preprocess_demand_collection!(
     # Allow the inflow from the primary network to be as large as required
     # (will be restricted when optimizing for the actual allocation)
     for link in p_independent.allocation.primary_network_connections[subnetwork_id]
-        JuMP.set_upper_bound(flow[link], MAX_ABS_FLOW / scaling.flow)
+        JuMP.set_upper_bound(
+            flow[link],
+            flow_capacity_upper_bound(link, p_independent) / scaling.flow,
+        )
         JuMP.set_lower_bound(flow[link], 0)
     end
 
@@ -360,14 +363,19 @@ function allocate_flows_to_subnetwork(
     primary_problem = primary_network.problem
 
     for secondary_network in get_secondary_networks(allocation_models)
-        (; problem, subnetwork_id) = secondary_network
-
-        for link in primary_network_connections[subnetwork_id]
-            primary_flow = primary_problem[:flow]
-            primary_flow_value = JuMP.value(primary_flow[link])
-            secondary_flow = problem[:flow]
-            JuMP.set_upper_bound(secondary_flow[link], primary_flow_value)
-            JuMP.set_lower_bound(secondary_flow[link], primary_flow_value)
+        for link in primary_network_connections[secondary_network.subnetwork_id]
+            allocated_flow_value =
+                JuMP.value(primary_problem[:flow][link]) * primary_network.scaling.flow /
+                secondary_network.scaling.flow
+
+            JuMP.set_upper_bound(
+                secondary_network.problem[:flow][link],
+                allocated_flow_value;
+            )
+            JuMP.set_lower_bound(
+                secondary_network.problem[:flow][link],
+                allocated_flow_value,
+            )
         end
     end
 end
@@ -444,40 +452,41 @@ function set_secondary_network_demands!(
     # Retrieve variable and constraint collections from the JuMP problem
     average_flow_unit_error = problem[:average_flow_unit_error]
     average_flow_unit_error_constraint = problem[:average_flow_unit_error_constraint]
+    for (demand_priority_idx, demand_priority) in enumerate(demand_priorities_all)
+        if !secondary_model.has_demand_priority[demand_priority_idx]
+            continue
+        end
+        # Objective metadata corresponding to this demand priority
+        (; expression_first) =
+            get_objective_data_of_demand_priority(objectives, demand_priority)
 
-    for link in keys(secondary_model.secondary_network_demand)
-        for (demand_priority_idx, demand_priority) in enumerate(demand_priorities_all)
-            if !secondary_model.has_demand_priority[demand_priority_idx]
-                continue
-            end
-            # Objective metadata corresponding to this demand priority
-            (; expression_first) =
-                get_objective_data_of_demand_priority(objectives, demand_priority)
-
-            demand = secondary_model.secondary_network_demand[link][demand_priority_idx]
+        for link in keys(secondary_model.secondary_network_demand)
+            d =
+                secondary_model.secondary_network_demand[link][demand_priority_idx] *
+                secondary_model.scaling.flow / primary_model.scaling.flow
 
-            # Demand is upper bound of what is allocated
-            JuMP.set_upper_bound(node_allocated[link, demand_priority], demand)
+            # Demand is upper bound of what can be allocated
+            JuMP.set_upper_bound(node_allocated[link, demand_priority], d)
 
             # Set demand in constraint for error term in first objective
             c = node_relative_error_constraint[link, demand_priority]
             error_term_first = node_error[link, demand_priority, :first]
-            JuMP.set_normalized_coefficient(c, error_term_first, demand)
-            JuMP.set_normalized_rhs(c, demand)
+            JuMP.set_normalized_coefficient(c, error_term_first, d)
+            JuMP.set_normalized_rhs(c, d)
 
             # Set demand in first objective expression
-            expression_first.terms[error_term_first] = demand
+            expression_first.terms[error_term_first] = d
 
             # Set demand in definition of average relative flow unit error
             JuMP.set_normalized_coefficient(
                 average_flow_unit_error_constraint[demand_priority],
                 error_term_first,
-                -demand,
+                -d,
             )
             add_to_coefficient!(
                 average_flow_unit_error_constraint[demand_priority],
                 average_flow_unit_error[demand_priority],
-                demand,
+                d,
             )
         end
     end
@@ -700,7 +709,12 @@ function optimize_multi_objective!(
         for link in primary_network_connections
             if type == AllocationObjectiveType.demand_flow ||
                type == AllocationObjectiveType.demand_storage
-                demand = JuMP.value(problem[:flow][link])
+                previous_demand = 0
+                if demand_priority_idx > 1
+                    previous_demand =
+                        secondary_model.secondary_network_demand[link][demand_priority_idx - 1]
+                end
+                demand = JuMP.value(problem[:flow][link]) - previous_demand
                 secondary_model.secondary_network_demand[link][demand_priority_idx] = demand
             end
         end
@@ -746,7 +760,6 @@ function optimize!(allocation_model::AllocationModel, model)::Nothing
                 "Allocation optimization for subnetwork $subnetwork_id at t = $t s is infeasible",
             )
         end
-    elseif termination_status != JuMP.OPTIMAL
     end
 
     return nothing
@@ -1039,7 +1052,8 @@ function update_allocation!(model)::Nothing
         set_simulation_data!(primary_network, integrator)
 
         reset_demand_coefficients(primary_network)
-        for secondary_network in get_secondary_networks(allocation_models)
+        for secondary_network in
+            sort(get_secondary_networks(allocation_models); by = x -> x.subnetwork_id)
             delete_temporary_constraints!(secondary_network)
             preprocess_demand_collection!(secondary_network, p_independent)
             optimize_multi_objective!(
@@ -1062,16 +1076,15 @@ function update_allocation!(model)::Nothing
         delete_temporary_constraints!(allocation_model)
         optimize!(allocation_model, model)
         parse_allocations!(integrator, allocation_model)
-
-        save_flows!(integrator, allocation_model, AllocationOptimizationType.allocate)
-        apply_control_from_allocation!(pump, allocation_model, integrator)
-        apply_control_from_allocation!(outlet, allocation_model, integrator)
-
         # allocate flows optimized from the primary network to the secondary networks
         if is_primary_network(allocation_model.subnetwork_id)
             allocate_flows_to_subnetwork(allocation_models, primary_network_connections)
         end
 
+        save_flows!(integrator, allocation_model, AllocationOptimizationType.allocate)
+        apply_control_from_allocation!(pump, allocation_model, integrator)
+        apply_control_from_allocation!(outlet, allocation_model, integrator)
+
         # Reset cumulative data
         reset_cumulative!(allocation_model)
     end

core/test/allocation_physics_test.jl

@@ -166,8 +166,7 @@ end
     @test ispath(toml_path)
 
     config = Ribasim.Config(toml_path; experimental_allocation = true)
-    model = Ribasim.Model(config)
-    Ribasim.solve!(model)
+    model = Ribasim.run(toml_path)
     allocation_flow_table = DataFrame(Ribasim.allocation_flow_data(model))
     basin_data = DataFrame(Ribasim.basin_data(model))
 
@@ -197,26 +196,42 @@ end
     using Ribasim
     using DataFrames: DataFrame
 
-    toml_path = normpath(
+    toml_path_1 = normpath(
         @__DIR__,
         "../../generated_testmodels/medium_primary_secondary_network/ribasim.toml",
     )
-    @test ispath(toml_path)
-
-    config = Ribasim.Config(toml_path; experimental_allocation = true)
-    model = Ribasim.run(toml_path)
-    allocation_flow_table = DataFrame(Ribasim.allocation_flow_data(model))
-    basin_data = DataFrame(Ribasim.basin_data(model))
-
-    link_outlet_3a = filter(:link_id => ==(23), allocation_flow_table)
-    link_outlet_3b = filter(:link_id => ==(24), allocation_flow_table)
-
-    flow_userdemand_primnet = filter(:link_id => ==(12), allocation_flow_table)
-    flow_userdemand_subnet_2 = filter(:link_id => ==(20), allocation_flow_table)
-    flow_userdemand_subnet_3 = filter(:link_id => ==(25), allocation_flow_table)
-
-    # Assert all 3 demands are met:
-    @test all(flow_userdemand_primnet.flow_rate .≈ 0.05)
-    @test all(flow_userdemand_subnet_2.flow_rate .≈ 0.1)
-    @test all(flow_userdemand_subnet_3.flow_rate .≈ 0.1)
+    toml_path_2 = normpath(
+        @__DIR__,
+        "../../generated_testmodels/medium_primary_secondary_network_verification/ribasim.toml",
+    )
+    model_1 = Ribasim.run(toml_path_1)
+    model_2 = Ribasim.run(toml_path_2)
+
+    flow_results_multiple_subnetwork = DataFrame(Ribasim.allocation_flow_data(model_1))
+    flow_results_single_subnetwork = DataFrame(Ribasim.allocation_flow_data(model_2))
+
+    #TODO: sometimes the model does not converge the first time step, rerun up to 5 times
+    rerun = 0
+    while length(filter(:link_id => ==(2), flow_results_multiple_subnetwork).flow_rate) !=
+          19 && rerun < 5
+        global rerun += 1
+        global model_2 = Ribasim.run(toml_path_2)
+        global flow_results_single_subnetwork =
+            DataFrame(Ribasim.allocation_flow_data(model_2))
+    end
+
+    # Assert that the flows over all links are the same
+    for link_id in unique(flow_results_multiple_subnetwork.link_id)
+        multiple_subs =
+            filter(:link_id => ==(link_id), flow_results_multiple_subnetwork).flow_rate
+        single_sub =
+            filter(:link_id => ==(link_id), flow_results_single_subnetwork).flow_rate
+        if !all(isapprox.(multiple_subs, single_sub; atol = 1e-8))
+            println(
+                "The flows over link $link_id differ by ",
+                maximum(single_sub .- multiple_subs),
+            )
+            @test false
+        end
+    end
 end

core/test/allocation_test.jl

@@ -623,21 +623,31 @@ end
         )
     end
 
-    @test_broken logger.logs[5].level == Error
-    @test_broken logger.logs[5].message == "Set of incompatible constraints found"
-    # @test sort(name.(keys(logger.logs[5].kwargs[:constraint_violations]))) == [
-    #     "linear_resistance_constraint[LinearResistance #2]",
-    #     "volume_conservation[Basin #1]",
-    # ]
-
-    # @test ispath(
-    #     @__DIR__,
-    #     "../../generated_testmodels/invalid_infeasible/results/allocation_analysis_infeasibility.log",
-    # )
-    # @test ispath(
-    #     @__DIR__,
-    #     "../../generated_testmodels/invalid_infeasible/results/allocation_analysis_scaling.log",
-    # )
+    # We need to check this way if this log message exists, because different hardware/OS can log in different order
+    found_message =
+        any(log -> log.message == "Set of incompatible constraints found", logger.logs)
+    @test found_message
+
+    constraint_names = [
+        "linear_resistance_constraint[LinearResistance #2]",
+        "volume_conservation[Basin #1]",
+    ]
+    found_message = any(
+        log ->
+            haskey(log.kwargs, :constraint_violations) &&
+                sort(name.(keys(log.kwargs[:constraint_violations]))) == constraint_names,
+        logger.logs,
+    )
+    @test found_message
+
+    @test ispath(
+        @__DIR__,
+        "../../generated_testmodels/invalid_infeasible/results/allocation_analysis_infeasibility.log",
+    )
+    @test ispath(
+        @__DIR__,
+        "../../generated_testmodels/invalid_infeasible/results/allocation_analysis_scaling.log",
+    )
 end
 
 @testitem "drain surplus" begin