Skip to content

Commit

Permalink
Begin to address input handling for lanl-ansi#99. Rename `head_curve_…
Browse files Browse the repository at this point in the history
…form` to `pump_type`.
  • Loading branch information
tasseff authored and hskkanth committed Feb 15, 2023
1 parent 0626c2e commit 4825842
Show file tree
Hide file tree
Showing 8 changed files with 75 additions and 62 deletions.
1 change: 1 addition & 0 deletions src/core/constants.jl
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,7 @@ JSON.lower(x::FLOW_DIRECTION) = Int(x)
PUMP_BEST_EFFICIENCY_POINT = 1 # Head gain takes a quadratic best efficiency form.
PUMP_EPANET = 2 # Head gain takes the form used by EPANET.
PUMP_LINEAR_POWER = 3 # Power is modeled linearly and head gain quadratic.
PUMP_CONSTANT_POWER = 4 # Power is constant when the pump is active.
end

"Ensures that JSON serialization of `PUMP` returns an integer."
Expand Down
8 changes: 4 additions & 4 deletions src/core/constraint_template.jl
Original file line number Diff line number Diff line change
Expand Up @@ -770,7 +770,7 @@ function constraint_on_off_pump_head_gain(
flow_transform(_FLOW_MIN),
)

if ref(wm, nw, :pump, a)["head_curve_form"] == PUMP_EPANET &&
if ref(wm, nw, :pump, a)["pump_type"] == PUMP_EPANET &&
isa(wm, AbstractNonlinearModel)
message = "PUMP_EPANET head curves are not currently supported for nonlinear models."
Memento.error(_LOGGER, message)
Expand Down Expand Up @@ -811,9 +811,9 @@ function constraint_on_off_pump_power(
_initialize_con_dict(wm, :on_off_pump_power, nw = nw, is_array = true)
con(wm, nw, :on_off_pump_power)[a] = Array{JuMP.ConstraintRef}([])

if ref(wm, nw, :pump, a)["head_curve_form"] in [PUMP_QUADRATIC, PUMP_EPANET]
if ref(wm, nw, :pump, a)["pump_type"] in [PUMP_QUADRATIC, PUMP_EPANET]
constraint_on_off_pump_power(wm, nw, a, q_min_forward)
elseif ref(wm, nw, :pump, a)["head_curve_form"] == PUMP_BEST_EFFICIENCY_POINT
elseif ref(wm, nw, :pump, a)["pump_type"] == PUMP_BEST_EFFICIENCY_POINT
wm_data = get_wm_data(wm.data)
density = _calc_scaled_density(wm_data)
gravity = _calc_scaled_gravity(wm_data)
Expand All @@ -825,7 +825,7 @@ function constraint_on_off_pump_power(
gravity,
q_min_forward,
)
elseif ref(wm, nw, :pump, a)["head_curve_form"] == PUMP_LINEAR_POWER
elseif ref(wm, nw, :pump, a)["pump_type"] == PUMP_LINEAR_POWER
# Ensure that the required keys for modeling pump power exist.
@assert haskey(ref(wm, nw, :pump, a), "power_fixed")
@assert haskey(ref(wm, nw, :pump, a), "power_per_unit_flow")
Expand Down
2 changes: 1 addition & 1 deletion src/core/data.jl
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,7 @@ end
function correct_enums!(data::Dict{String,<:Any})
correct_statuses!(data)
correct_flow_directions!(data)
correct_pump_head_curve_forms!(data)
correct_pump_types!(data)
end


Expand Down
56 changes: 28 additions & 28 deletions src/core/pump.jl
Original file line number Diff line number Diff line change
Expand Up @@ -56,24 +56,24 @@ function set_pump_flow_partition!(
end


function correct_pump_head_curve_forms!(data::Dict{String,<:Any})
apply_wm!(_correct_pump_head_curve_forms!, data; apply_to_subnetworks = true)
function correct_pump_types!(data::Dict{String,<:Any})
apply_wm!(_correct_pump_types!, data; apply_to_subnetworks = true)
end


function _correct_pump_head_curve_forms!(data::Dict{String,<:Any})
function _correct_pump_types!(data::Dict{String,<:Any})
components = values(get(data, "pump", Dict{String,Any}()))
_correct_pump_head_curve_form!.(components)
_correct_pump_type!.(components)
end


function _correct_pump_head_curve_form!(pump::Dict{String, <:Any})
head_curve_tmp = get(pump, "head_curve_form", PUMP_QUADRATIC)
function _correct_pump_type!(pump::Dict{String, <:Any})
head_curve_tmp = get(pump, "pump_type", PUMP_QUADRATIC)

if isa(head_curve_tmp, PUMP)
pump["head_curve_form"] = head_curve_tmp
pump["pump_type"] = head_curve_tmp
else
pump["head_curve_form"] = PUMP(head_curve_tmp)
pump["pump_type"] = PUMP(head_curve_tmp)
end
end

Expand Down Expand Up @@ -108,10 +108,10 @@ function _calc_pump_head_gain_max(pump::Dict{String, <:Any}, node_fr::Dict{Strin
# Calculate the flow at the maximum head gain, then return maximum head gain.
c = _calc_head_curve_coefficients(pump)

if pump["head_curve_form"] in [PUMP_QUADRATIC, PUMP_BEST_EFFICIENCY_POINT, PUMP_LINEAR_POWER]
if pump["pump_type"] in [PUMP_QUADRATIC, PUMP_BEST_EFFICIENCY_POINT, PUMP_LINEAR_POWER]
flow_at_max = -c[2] * inv(2.0 * c[1]) > 0.0 ? -c[2] * inv(2.0 * c[1]) : 0.0
return max(0.0, c[1] * flow_at_max^2 + c[2] * flow_at_max + c[3])
elseif pump["head_curve_form"] == PUMP_EPANET
elseif pump["pump_type"] == PUMP_EPANET
return max(0.0, c[1])
end
end
Expand All @@ -121,7 +121,7 @@ function _calc_pump_flow_max(pump::Dict{String,<:Any}, node_fr::Dict{String,Any}
# Get possible maximal flow values based on the head curve.
c = _calc_head_curve_coefficients(pump)

if pump["head_curve_form"] in [PUMP_QUADRATIC, PUMP_BEST_EFFICIENCY_POINT, PUMP_LINEAR_POWER]
if pump["pump_type"] in [PUMP_QUADRATIC, PUMP_BEST_EFFICIENCY_POINT, PUMP_LINEAR_POWER]
q_max_1 = (-c[2] + sqrt(c[2]^2 - 4.0*c[1]*c[3])) * inv(2.0 * c[1])
q_max_2 = (-c[2] - sqrt(c[2]^2 - 4.0*c[1]*c[3])) * inv(2.0 * c[1])

Expand All @@ -132,7 +132,7 @@ function _calc_pump_flow_max(pump::Dict{String,<:Any}, node_fr::Dict{String,Any}

# Get the minimal value of the above and the possible "flow_max" value.
return min(max(q_max_1, q_max_2), max(q_max_3, q_max_4), get(pump, "flow_max", Inf))
elseif pump["head_curve_form"] == PUMP_EPANET
elseif pump["pump_type"] == PUMP_EPANET
return min((-c[1] * inv(c[2]))^(inv(c[3])), get(pump, "flow_max", Inf))
end
end
Expand Down Expand Up @@ -160,61 +160,61 @@ end


function _calc_head_curve_coefficients(pump::Dict{String, <:Any})
if pump["head_curve_form"] in [PUMP_QUADRATIC, PUMP_LINEAR_POWER]
if pump["pump_type"] in [PUMP_QUADRATIC, PUMP_LINEAR_POWER]
return _calc_head_curve_coefficients_quadratic(pump)
elseif pump["head_curve_form"] == PUMP_BEST_EFFICIENCY_POINT
elseif pump["pump_type"] == PUMP_BEST_EFFICIENCY_POINT
return _calc_head_curve_coefficients_best_efficiency_point(pump)
elseif pump["head_curve_form"] == PUMP_EPANET
elseif pump["pump_type"] == PUMP_EPANET
return _calc_head_curve_coefficients_epanet(pump)
else
error("\"$(pump["head_curve_form"])\" is not a valid head curve formulation.")
error("\"$(pump["pump_type"])\" is not a valid head curve formulation.")
end
end


function _calc_head_curve_function(pump::Dict{String, <:Any})
if pump["head_curve_form"] in [PUMP_QUADRATIC, PUMP_LINEAR_POWER]
if pump["pump_type"] in [PUMP_QUADRATIC, PUMP_LINEAR_POWER]
coeff = _calc_head_curve_coefficients_quadratic(pump)
return x -> sum(coeff .* [x^2, x, 1.0])
elseif pump["head_curve_form"] == PUMP_BEST_EFFICIENCY_POINT
elseif pump["pump_type"] == PUMP_BEST_EFFICIENCY_POINT
coeff = _calc_head_curve_coefficients_best_efficiency_point(pump)
return x -> sum(coeff .* [x^2, x, 1.0])
elseif pump["head_curve_form"] == PUMP_EPANET
elseif pump["pump_type"] == PUMP_EPANET
coeff = _calc_head_curve_coefficients_epanet(pump)
return x -> coeff[1] + coeff[2] * x^coeff[3]
else
error("\"$(pump["head_curve_form"])\" is not a valid head curve formulation.")
error("\"$(pump["pump_type"])\" is not a valid head curve formulation.")
end
end


function _calc_head_curve_function(pump::Dict{String, <:Any}, z::JuMP.VariableRef)
if pump["head_curve_form"] in [PUMP_QUADRATIC, PUMP_LINEAR_POWER]
if pump["pump_type"] in [PUMP_QUADRATIC, PUMP_LINEAR_POWER]
coeff = _calc_head_curve_coefficients_quadratic(pump)
return x -> sum(coeff .* [x^2, x, z])
elseif pump["head_curve_form"] == PUMP_BEST_EFFICIENCY_POINT
elseif pump["pump_type"] == PUMP_BEST_EFFICIENCY_POINT
coeff = _calc_head_curve_coefficients_best_efficiency_point(pump)
return x -> sum(coeff .* [x^2, x, z])
elseif pump["head_curve_form"] == PUMP_EPANET
elseif pump["pump_type"] == PUMP_EPANET
coeff = _calc_head_curve_coefficients_epanet(pump)
return x -> coeff[1] * z + coeff[2] * x^coeff[3]
else
error("\"$(pump["head_curve_form"])\" is not a valid head curve formulation.")
error("\"$(pump["pump_type"])\" is not a valid head curve formulation.")
end
end

function _calc_head_curve_derivative(pump::Dict{String, <:Any})
if pump["head_curve_form"] in [PUMP_QUADRATIC, PUMP_LINEAR_POWER]
if pump["pump_type"] in [PUMP_QUADRATIC, PUMP_LINEAR_POWER]
coeff = _calc_head_curve_coefficients_quadratic(pump)
return x -> sum(coeff .* [2.0 * x, 1.0, 0.0])
elseif pump["head_curve_form"] == PUMP_BEST_EFFICIENCY_POINT
elseif pump["pump_type"] == PUMP_BEST_EFFICIENCY_POINT
coeff = _calc_head_curve_coefficients_best_efficiency_point(pump)
return x -> sum(coeff .* [2.0 * x, 1.0, 0.0])
elseif pump["head_curve_form"] == PUMP_EPANET
elseif pump["pump_type"] == PUMP_EPANET
coeff = _calc_head_curve_coefficients_epanet(pump)
return x -> coeff[2] * coeff[3] * x^(coeff[3] - 1.0)
else
error("\"$(pump["head_curve_form"])\" is not a valid head curve formulation.")
error("\"$(pump["pump_type"])\" is not a valid head curve formulation.")
end
end

Expand Down
58 changes: 35 additions & 23 deletions src/io/epanet.jl
Original file line number Diff line number Diff line change
Expand Up @@ -922,33 +922,45 @@ function _read_pump!(data::Dict{String,<:Any})
# Loop over remaining entries and store remaining properties.
for i in range(4; stop = length(current), step = 2)
if uppercase(current[i]) == "POWER"
Memento.error(_LOGGER, "Constant power pumps are not supported.")
elseif uppercase(current[i]) != "HEAD"
Memento.error(_LOGGER, "Pump keyword in INP file not recognized.")
end
# Memento.error(_LOGGER, "Constant power pumps are not supported.")
pump["pump_type"] = PUMP_CONSTANT_POWER

# Obtain and scale head-versus-flow pump curve.
flow = first.(data["curve"][current[i+1]]) # Flow.
head = last.(data["curve"][current[i+1]]) # Head.
# Parse the head of the reservoir (in meters).
if data["flow_units"] == "LPS" || data["flow_units"] == "CMH"
# Conver power from kilowatts to watts.
pump["power"] = 1.0e3 * parse(Float64, current[i+1])
elseif data["flow_units"] == "GPM" # If gallons per minute...
# Convert power from horsepower to watts.
pump["power"] = 745.7 * parse(Float64, current[i+1])
else
Memento.error(_LOGGER, "Could not find a valid \"units\" option type.")
end
elseif uppercase(current[i]) == "HEAD"
# Obtain and scale head-versus-flow pump curve.
flow = first.(data["curve"][current[i+1]]) # Flow.
head = last.(data["curve"][current[i+1]]) # Head.

if data["flow_units"] == "LPS" # If liters per second...
# Convert from liters per second to cubic meters per second.
flow *= 1.0e-3
elseif data["flow_units"] == "CMH" # If cubic meters per hour...
# Convert from cubic meters per hour to cubic meters per second.
flow *= inv(3600.0)
elseif data["flow_units"] == "GPM" # If gallons per minute...
# Convert from gallons per minute to cubic meters per second.
flow *= 6.30902e-5
if data["flow_units"] == "LPS" # If liters per second...
# Convert from liters per second to cubic meters per second.
flow *= 1.0e-3
elseif data["flow_units"] == "CMH" # If cubic meters per hour...
# Convert from cubic meters per hour to cubic meters per second.
flow *= inv(3600.0)
elseif data["flow_units"] == "GPM" # If gallons per minute...
# Convert from gallons per minute to cubic meters per second.
flow *= 6.30902e-5

# Convert head from feet to meters.
head *= 0.3048
else
Memento.error(_LOGGER, "Could not find a valid \"units\" option type.")
end
# Convert head from feet to meters.
head *= 0.3048
else
Memento.error(_LOGGER, "Could not find a valid \"units\" option type.")
end

# Curve of head (meters) versus flow (cubic meters per second).
pump["head_curve"] = Array([(flow[j], head[j]) for j = 1:length(flow)])
# Curve of head (meters) versus flow (cubic meters per second).
pump["head_curve"] = Array([(flow[j], head[j]) for j = 1:length(flow)])
elseif uppercase(current[i]) != "HEAD"
Memento.error(_LOGGER, "Pump keyword in INP file not recognized.")
end
end

# Flow is always in the positive direction for pumps.
Expand Down
2 changes: 1 addition & 1 deletion test/io.jl
Original file line number Diff line number Diff line change
Expand Up @@ -37,6 +37,6 @@
network = WaterModels.parse_file("../test/data/epanet/snapshot/pump-hw-lps.inp")
data = JSON.parse(JSON.json(network)) # Convert data to JSON dictionary.
@test data["pump"]["1"]["flow_direction"] == 1
@test data["pump"]["1"]["head_curve_form"] == 0
@test data["pump"]["1"]["pump_type"] == 0
end
end
4 changes: 2 additions & 2 deletions test/owf.jl
Original file line number Diff line number Diff line change
Expand Up @@ -40,7 +40,7 @@ for formulation in [NCWaterModel, NCDWaterModel, CRDWaterModel, LAWaterModel, LR
@testset "Optimal Water Flow Problems (Single Network, EPANET Pump Formulation): $(formulation)" begin
network = WaterModels.parse_file("../test/data/epanet/snapshot/pump-hw-lps.inp")

map(x -> x["head_curve_form"] = PUMP_EPANET, values(network["pump"]))
map(x -> x["pump_type"] = PUMP_EPANET, values(network["pump"]))
WaterModels.recompute_bounds!(network) # Recompute component bounds after the above changes.
set_flow_partitions_si!(network, 10.0, 1.0e-4)

Expand Down Expand Up @@ -90,4 +90,4 @@ end
result = WaterModels.solve_mn_owf(network_mn, LRDWaterModel, milp_solver)
result = WaterModels.run_mn_owf(network_mn, LRDWaterModel, milp_solver)
@test result["termination_status"] == OPTIMAL
end
end
6 changes: 3 additions & 3 deletions test/pump.jl
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
@testset "src/core/pump.jl" begin
for head_curve_form in [PUMP_QUADRATIC, PUMP_BEST_EFFICIENCY_POINT, PUMP_EPANET, PUMP_LINEAR_POWER]
@testset "pump with efficiency curve: $(head_curve_form)" begin
for pump_type in [PUMP_QUADRATIC, PUMP_BEST_EFFICIENCY_POINT, PUMP_EPANET, PUMP_LINEAR_POWER]
@testset "pump with efficiency curve: $(pump_type)" begin
# Read in the initial network data and ensure an efficiency curve exists.
network_path = "../test/data/epanet/snapshot/pump-hw-lps-eff-curve.inp"
network = WaterModels.parse_file(network_path)
Expand All @@ -21,7 +21,7 @@
map(x -> x["power_per_unit_flow"] = expected_val, values(network["pump"]))

# Change the head curve form and recompute network bounds.
map(x -> x["head_curve_form"] = head_curve_form, values(network["pump"]))
map(x -> x["pump_type"] = pump_type, values(network["pump"]))
WaterModels.recompute_bounds!(network)

# Set up and solve an optimization problem with the pump data.
Expand Down

0 comments on commit 4825842

Please sign in to comment.