Skip to content

Commit

Permalink
Comments adressed
Browse files Browse the repository at this point in the history
  • Loading branch information
SouthEndMusic committed Apr 11, 2024
1 parent ff1cf3a commit ab91921
Show file tree
Hide file tree
Showing 3 changed files with 26 additions and 65 deletions.
6 changes: 3 additions & 3 deletions core/src/allocation_optim.jl
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
@enumx OptimizationType internal_sources collect_demands allocate

"""
Add an objective term (1 - flow/demand)^2. If the absolute
Add an objective term `demand * (1 - flow/demand)^2`. If the absolute
value of the demand is very small, this would lead to huge coefficients,
so in that case a term of the form (flow - demand)^2 is used.
"""
Expand All @@ -11,12 +11,12 @@ function add_objective_term!(
F::JuMP.VariableRef,
)::Nothing
if abs(demand) < 1e-5
# Error term (F - d)^2
# Error term (F - d)^2 = F² - 2dF + d²
JuMP.add_to_expression!(ex, 1.0, F, F)
JuMP.add_to_expression!(ex, -2.0 * demand, F)
JuMP.add_to_expression!(ex, demand^2)
else
# Error term d*(1 - F/d)^2
# Error term d*(1 - F/d)^2 = F²/d - 2F + d
JuMP.add_to_expression!(ex, 1.0 / demand, F, F)
JuMP.add_to_expression!(ex, -2.0, F)
JuMP.add_to_expression!(ex, demand)
Expand Down
83 changes: 22 additions & 61 deletions core/test/allocation_test.jl
Original file line number Diff line number Diff line change
Expand Up @@ -223,21 +223,11 @@ end

# See the difference between these values here and in
# "subnetworks_with_sources"
@test all(
isapprox.(
subnetwork_demands[(NodeID(:Basin, 2), NodeID(:Pump, 11))],
[4.0, 4.0, 0.0],
atol = 1e-4,
),
)
@test subnetwork_demands[(NodeID(:Basin, 2), NodeID(:Pump, 11))] [4.0, 4.0, 0.0] atol =
1e-4
@test subnetwork_demands[(NodeID(:Basin, 6), NodeID(:Pump, 24))] [0.004, 0.0, 0.0]
@test all(
isapprox.(
subnetwork_demands[(NodeID(:Basin, 10), NodeID(:Pump, 38))][1:2],
[0.001, 0.002],
rtol = 1e-4,
),
)
@test subnetwork_demands[(NodeID(:Basin, 10), NodeID(:Pump, 38))][1:2] [0.001, 0.002] rtol =
1e-4

# Solving for the main network, containing subnetworks as UserDemands
allocation_model = allocation_models[1]
Expand All @@ -260,25 +250,13 @@ end
u = ComponentVector(; storage = zeros(length(p.basin.node_id)))
Ribasim.update_allocation!((; p, t, u))

@test all(
isapprox.(
subnetwork_allocateds[NodeID(:Basin, 2), NodeID(:Pump, 11)],
[4.0, 0.4947, 0.0],
atol = 1e-4,
),
)
@test isapprox(
subnetwork_allocateds[NodeID(:Basin, 6), NodeID(:Pump, 24)],
[0.004, 0.0, 0.0],
rtol = 1e-3,
)
@test all(
isapprox.(
subnetwork_allocateds[NodeID(:Basin, 10), NodeID(:Pump, 38)],
[0.001, 2.473e-4, 0.0],
rtol = 1e-3,
),
)
@test subnetwork_allocateds[NodeID(:Basin, 2), NodeID(:Pump, 11)] [4.0, 0.4947, 0.0] atol =
1e-4
@test subnetwork_allocateds[NodeID(:Basin, 6), NodeID(:Pump, 24)] [0.004, 0.0, 0.0] rtol =
1e-3

@test subnetwork_allocateds[NodeID(:Basin, 10), NodeID(:Pump, 38)]
[0.001, 2.473e-4, 0.0] rtol = 1e-3

# Test for existence of edges in allocation flow record
allocation_flow = DataFrame(record_flow)
Expand Down Expand Up @@ -322,21 +300,15 @@ end
# Collecting demands
u = ComponentVector(; storage = zeros(length(basin.node_id)))
for (i, allocation_model) in enumerate(allocation_models[2:end])
@show i
Ribasim.allocate!(p, allocation_model, t, u, OptimizationType.internal_sources)
Ribasim.allocate!(p, allocation_model, t, u, OptimizationType.collect_demands)
end

# See the difference between these values here and in
# "allocation with main network optimization problem", internal sources
# lower the subnetwork demands
@test all(
isapprox.(
subnetwork_demands[(NodeID(:Basin, 2), NodeID(:Pump, 11))],
[3.1, 4.0, 0.0],
atol = 1e-4,
),
)
@test subnetwork_demands[(NodeID(:Basin, 2), NodeID(:Pump, 11))] [3.1, 4.0, 0.0] atol =
1e-4
@test subnetwork_demands[(NodeID(:Basin, 6), NodeID(:Pump, 24))] [0.004, 0.0, 0.0]
@test subnetwork_demands[(NodeID(:Basin, 10), NodeID(:Pump, 38))][1:2] [0.001, 0.001]
end
Expand Down Expand Up @@ -471,7 +443,7 @@ end
)
objective = JuMP.objective_function(problem)
# Reduced demand
@test isapprox(flow_demand.demand[1], flow_demand.demand_itp[1](t) - 0.001, rtol = 1e-3)
@test flow_demand.demand[1] flow_demand.demand_itp[1](t) - 0.001 rtol = 1e-3
@test JuMP.normalized_rhs(constraint_flow_out) == Inf

## Priority 2
Expand All @@ -486,7 +458,7 @@ end
# No demand left
@test flow_demand.demand[1] < 1e-10
# Allocated
@test isapprox(JuMP.value(only(F_flow_buffer_in)), 0.001, rtol = 1e-3)
@test JuMP.value(only(F_flow_buffer_in)) 0.001 rtol = 1e-3
@test JuMP.normalized_rhs(constraint_flow_out) == 0.0

## Priority 3
Expand All @@ -502,17 +474,10 @@ end
# The flow from the source is used up in previous priorities
@test JuMP.value(F[(NodeID(NodeType.LevelBoundary, 1), node_id_with_flow_demand)]) == 0
# So flow from the flow buffer is used for UserDemand #4
@test isapprox(
JuMP.value(F_flow_buffer_out[node_id_with_flow_demand]),
0.001,
rtol = 1e-3,
)
@test JuMP.value(F_flow_buffer_out[node_id_with_flow_demand]) 0.001 rtol = 1e-3

# Flow taken from buffer
@test isapprox(
JuMP.value(only(F_flow_buffer_out)),
user_demand.demand_itp[1][3](t),
rtol = 1e-3,
)
@test JuMP.value(only(F_flow_buffer_out)) user_demand.demand_itp[1][3](t) rtol = 1e-3
# No flow coming from level boundary
@test JuMP.value(F[(only(level_boundary.node_id), node_id_with_flow_demand)]) == 0

Expand All @@ -527,12 +492,9 @@ end
)
# Get demand from buffers
d = user_demand.demand_itp[3][4](t)
@test isapprox(
JuMP.value(F[(NodeID(NodeType.UserDemand, 4), NodeID(NodeType.Basin, 7))]) +
JuMP.value(F[(NodeID(NodeType.UserDemand, 6), NodeID(NodeType.Basin, 7))]),
d,
rtol = 1e-3,
)
@test JuMP.value(F[(NodeID(NodeType.UserDemand, 4), NodeID(NodeType.Basin, 7))]) +
JuMP.value(F[(NodeID(NodeType.UserDemand, 6), NodeID(NodeType.Basin, 7))]) d rtol =
1e-3
end

@testitem "flow_demand_with_max_flow_rate" begin
Expand Down Expand Up @@ -572,8 +534,7 @@ end
fractions = Vector{Float64}[]

for id in user_demand.node_id
data_allocation_id =
filter(:node_id => node_id -> node_id == id.value, data_allocation)
data_allocation_id = filter(:node_id => ==(id.value), data_allocation)
frac = data_allocation_id.allocated ./ data_allocation_id.demand
push!(fractions, frac)
end
Expand Down
2 changes: 1 addition & 1 deletion python/ribasim_testmodels/ribasim_testmodels/allocation.py
Original file line number Diff line number Diff line change
Expand Up @@ -1034,7 +1034,7 @@ def linear_resistance_demand_model():
def fair_distribution_model():
model = Model(
starttime="2020-01-01 00:00:00",
endtime="2021-01-01 00:00:00",
endtime="2020-01-07 00:00:00",
crs="EPSG:28992",
allocation=Allocation(use_allocation=True),
)
Expand Down

0 comments on commit ab91921

Please sign in to comment.