diff --git a/src/WaterModels.jl b/src/WaterModels.jl index 360043c5..07787850 100644 --- a/src/WaterModels.jl +++ b/src/WaterModels.jl @@ -86,6 +86,7 @@ include("form/crd.jl") include("form/la.jl") include("form/pwlrd.jl") include("form/lrd.jl") +include("form/scip_constraints_variables.jl") include("prob/wf.jl") include("prob/owf.jl") @@ -98,6 +99,7 @@ include("util/variable_index.jl") include("util/bound_problem.jl") include("util/pairwise_cuts.jl") include("util/obbt.jl") +include("util/fixing_variables.jl") # Deprecated functions. include("deprecated.jl") diff --git a/src/core/data.jl b/src/core/data.jl index dddb241f..8e5689ec 100644 --- a/src/core/data.jl +++ b/src/core/data.jl @@ -72,17 +72,19 @@ function _check_connectivity(data::Dict{String,<:Any}) end for comp_type in _LINK_COMPONENTS - for (i, comp) in data[comp_type] - if !(comp["node_fr"] in node_ids) - error_message = "From node $(comp["node_fr"]) in " - error_message *= "$(replace(comp_type, "_" => " ")) $(i) is not defined." - Memento.error(_LOGGER, error_message) - end + if(haskey(data,comp_type)) + for (i, comp) in data[comp_type] + if !(comp["node_fr"] in node_ids) + error_message = "From node $(comp["node_fr"]) in " + error_message *= "$(replace(comp_type, "_" => " ")) $(i) is not defined." + Memento.error(_LOGGER, error_message) + end - if !(comp["node_to"] in node_ids) - error_message = "To node $(comp["node_to"]) in " - error_message *= "$(replace(comp_type, "_" => " ")) $(i) is not defined." - Memento.error(_LOGGER, error_message) + if !(comp["node_to"] in node_ids) + error_message = "To node $(comp["node_to"]) in " + error_message *= "$(replace(comp_type, "_" => " ")) $(i) is not defined." + Memento.error(_LOGGER, error_message) + end end end end @@ -112,17 +114,19 @@ function _check_status(data::Dict{String,<:Any}) end for comp_type in _LINK_COMPONENTS - for (i, comp) in data[comp_type] - if comp["status"] != STATUS_INACTIVE && !(comp["node_fr"] in active_node_ids) - warning_message = "Active $(comp_type) $(i) is connected to inactive " - warning_message *= "from node $(comp["node_fr"])." - Memento.warn(_LOGGER, warning_message) - end + if(haskey(data,comp_type)) + for (i, comp) in data[comp_type] + if comp["status"] != STATUS_INACTIVE && !(comp["node_fr"] in active_node_ids) + warning_message = "Active $(comp_type) $(i) is connected to inactive " + warning_message *= "from node $(comp["node_fr"])." + Memento.warn(_LOGGER, warning_message) + end - if comp["status"] != STATUS_INACTIVE && !(comp["node_to"] in active_node_ids) - warning_message = "Active $(comp_type) $(i) is connected to inactive " - warning_message *= "to node $(comp["node_to"])." - Memento.warn(_LOGGER, warning_message) + if comp["status"] != STATUS_INACTIVE && !(comp["node_to"] in active_node_ids) + warning_message = "Active $(comp_type) $(i) is connected to inactive " + warning_message *= "to node $(comp["node_to"])." + Memento.warn(_LOGGER, warning_message) + end end end end @@ -512,6 +516,14 @@ function _calc_head_max(data::Dict{String, <:Any}) head_max = max(head_max, node_to["elevation"] + head_gain) end + for (i, pump) in wm_data["ne_pump"] + # Consider possible pump head gains in computation of head_max. + node_fr = wm_data["node"][string(pump["node_fr"])] + node_to = wm_data["node"][string(pump["node_to"])] + head_gain = _calc_pump_head_gain_max(pump, node_fr, node_to) + head_max = max(head_max, node_to["elevation"] + head_gain) + end + for (i, regulator) in wm_data["regulator"] # Consider possible downstream regulator heads in computation of head_max. p_setting, node_to_index = regulator["setting"], string(regulator["node_to"]) @@ -1023,7 +1035,7 @@ end function _apply_pipe_unit_transform!(data::Dict{String,<:Any}, transform_length::Function, head_loss::String) wm_data = get_wm_data(data) - + if !haskey(wm_data, "pipe") return end @@ -1139,6 +1151,78 @@ function _apply_pump_unit_transform!( end end +function _apply_ne_pump_unit_transform!( + data::Dict{String,<:Any}, transform_mass::Function, transform_flow::Function, + transform_length::Function, transform_time::Function) + wm_data = get_wm_data(data) + + if !haskey(wm_data, "ne_pump") + return + end + + power_scalar = transform_mass(1.0) * transform_length(1.0)^2 / transform_time(1.0)^3 + energy_scalar = transform_mass(1.0) * transform_length(1.0)^2 / transform_time(1.0)^2 + + for (_, ne_pump) in wm_data["ne_pump"] + if haskey(ne_pump, "head_curve") + ne_pump["head_curve"] = [(transform_flow(x[1]), x[2]) for x in ne_pump["head_curve"]] + ne_pump["head_curve"] = [(x[1], transform_length(x[2])) for x in ne_pump["head_curve"]] + end + + if haskey(ne_pump, "efficiency_curve") + ne_pump["efficiency_curve"] = [(transform_flow(x[1]), x[2]) for x in ne_pump["efficiency_curve"]] + end + + if haskey(ne_pump, "energy_price") + ne_pump["energy_price"] /= energy_scalar + end + + if haskey(ne_pump, "E") + ne_pump["E"] *= energy_scalar + end + + if haskey(ne_pump, "power_fixed") + ne_pump["power_fixed"] *= power_scalar + end + + if haskey(ne_pump, "P") + ne_pump["P"] *= power_scalar + end + + if haskey(ne_pump, "c") + ne_pump["c"] *= energy_scalar + end + + if haskey(ne_pump, "min_inactive_time") + ne_pump["min_inactive_time"] = transform_time(ne_pump["min_inactive_time"]) + end + + if haskey(ne_pump, "min_active_time") + ne_pump["min_active_time"] = transform_time(ne_pump["min_active_time"]) + end + + if haskey(ne_pump, "power_per_unit_flow") + ne_pump["power_per_unit_flow"] *= power_scalar / transform_flow(1.0) + end + end + + if haskey(wm_data, "time_series") && haskey(wm_data["time_series"], "ne_pump") + for ne_pump in values(wm_data["time_series"]["ne_pump"]) + if haskey(ne_pump, "energy_price") + ne_pump["energy_price"] ./= energy_scalar + end + + if haskey(ne_pump, "power_fixed") + ne_pump["power_fixed"] .*= power_scalar + end + + if haskey(ne_pump, "power_per_unit_flow") + ne_pump["power_per_unit_flow"] .*= power_scalar / transform_flow(1.0) + end + end + end +end + function _apply_regulator_unit_transform!(data::Dict{String,<:Any}, transform_head::Function, transform_length::Function) wm_data = get_wm_data(data) @@ -1169,7 +1253,7 @@ function _calc_scaled_gravity(data::Dict{String, <:Any}) if wm_data["per_unit"] base_time = 1.0 / _calc_time_per_unit_transform(wm_data)(1.0) - base_length = 1.0 / _calc_length_per_unit_transform(wm_data)(1.0) + base_length = 1.0 / _calc_length_per_unit_transform(wm_data)(1.0) return _calc_scaled_gravity(base_length, base_time) else return _GRAVITY @@ -1368,6 +1452,8 @@ function _transform_nw!( _apply_des_pipe_unit_transform!(wm_nw_data, length_transform, head_loss) _apply_pump_unit_transform!(wm_nw_data, mass_transform, flow_transform, length_transform, time_transform) + _apply_ne_pump_unit_transform!(wm_nw_data, mass_transform, + flow_transform, length_transform, time_transform) _apply_regulator_unit_transform!(wm_nw_data, head_transform, length_transform) # Apply transformations to nodal components. diff --git a/src/core/function.jl b/src/core/function.jl index b281fd47..64fe5c1c 100644 --- a/src/core/function.jl +++ b/src/core/function.jl @@ -1,84 +1,27 @@ -################################################################################# -# This file defines the nonlinear head loss functions for water systems models. # -################################################################################# -function _f_alpha(alpha::Float64; convex::Bool = false) - if convex - return function (x::Float64) - return (x * x)^(0.5 + 0.5 * alpha) - end - else - return function (x::Float64) - return sign(x) * (x * x)^(0.5 + 0.5 * alpha) - end - end -end - +# This file defines the nonlinear head loss functions for water systems models. # +################################################################################# -function _df_alpha(alpha::Float64; convex::Bool = false) - if convex - return function (x::Float64) - return (1.0 + alpha) * sign(x) * (x * x)^(0.5 * alpha) - end - else - return function (x::Float64) - return (1.0 + alpha) * (x * x)^(0.5 * alpha) - end +function _get_alpha(wm::AbstractWaterModel) + head_loss_form = wm.ref[:it][wm_it_sym][:head_loss] + alphas = uppercase(head_loss_form) == "H-W" ? 1.852 : 2.0 + # alphas = [ref(wm, nw, :alpha) for nw in nw_ids(wm)] + if any(x -> x != alphas[1], alphas) + Memento.error(_LOGGER, "Head loss exponents are different across the multinetwork.") end + return alphas[1] end - -function _d2f_alpha(alpha::Float64; convex::Bool = false) - if convex - return function (x::Float64) - return alpha * (1.0 + alpha) * (x * x)^(0.5 * alpha - 0.5) - end - else - return function (x::Float64) - return x != 0.0 ? - sign(x) * alpha * (1.0 + alpha) * (x * x)^(0.5 * alpha - 0.5) : - prevfloat(Inf) - end - end +function head_loss(wm::Union{NCDWaterModel,CRDWaterModel}, x) + p = _get_alpha(wm) + # An expression equivalent to abspower(x, p) = abs(x)^p + # return JuMP.@expression(wm.model, ifelse(x > 0, x^p, (-x)^p)) + return JuMP.@expression(wm.model, (x^2)^(p / 2)) end - -function _get_alpha_min_1(wm::AbstractWaterModel) - return _get_exponent_from_head_loss_form(wm.ref[:it][wm_it_sym][:head_loss]) - 1.0 -end - - -function head_loss_args(wm::Union{NCDWaterModel,CRDWaterModel}) - alpha_m1 = _get_alpha_min_1(wm) - return ( - :head_loss, - 1, - _f_alpha(alpha_m1, convex = true), - _df_alpha(alpha_m1, convex = true), - _d2f_alpha(alpha_m1, convex = true), - ) -end - - -function head_loss_args(wm::NCWaterModel) - alpha_m1 = _get_alpha_min_1(wm) - - return ( - :head_loss, - 1, - _f_alpha(alpha_m1, convex = false), - _df_alpha(alpha_m1, convex = false), - _d2f_alpha(alpha_m1, convex = false), - ) -end - - -function _function_head_loss(wm::AbstractWaterModel) - # By default, head loss is not defined by nonlinear registered functions. -end - - -function _function_head_loss(wm::AbstractNonlinearModel) - JuMP.register(wm.model, head_loss_args(wm)...) +function head_loss(wm::AbstractNonlinearModel, x) + p = _get_alpha(wm) + # An expression equivalent to signpower(x, p) = sign(x) * abs(x)^p + return JuMP.@expression(wm.model, ifelse(x > 0, x^p, -(-x)^p)) end diff --git a/src/core/objective.jl b/src/core/objective.jl index 1766eaf3..44a52761 100644 --- a/src/core/objective.jl +++ b/src/core/objective.jl @@ -85,32 +85,25 @@ cost of building and operating all network expansion components is minimized. """ function objective_ne(wm::AbstractWaterModel)::JuMP.AffExpr # Get all network IDs in the multinetwork. - network_ids = sort(collect(nw_ids(wm))) + first_network_id = sort(collect(nw_ids(wm)))[1] - # Find the network IDs over which the objective will be defined. - if length(network_ids) > 1 - network_ids_flow = network_ids[1:end-1] - else - network_ids_flow = network_ids - end # Initialize the objective expression to zero. objective = JuMP.AffExpr(0.0) - for n in network_ids_flow - # Get the set of network expansion short pipes at time index `n`. - for (a, ne_short_pipe) in ref(wm, n, :ne_short_pipe) - # Add the cost of network expansion component `a` at time period `n`. - term = ne_short_pipe["construction_cost"] * var(wm, n, :z_ne_short_pipe, a) - JuMP.add_to_expression!(objective, term) - end +n = first_network_id + # Get the set of network expansion short pipes at time index `n`. + for (a, ne_short_pipe) in ref(wm, n, :ne_short_pipe) + # Add the cost of network expansion component `a` at time period `n`. + term = ne_short_pipe["construction_cost"] * var(wm, n, :z_ne_short_pipe, a) + JuMP.add_to_expression!(objective, term) + end - # Get the set of network expansion pumps at time index `n`. - for (a, ne_pump) in ref(wm, n, :ne_pump) - # Add the cost of network expansion component `a` at time period `n`. - term = ne_pump["construction_cost"] * var(wm, n, :z_ne_pump, a) - JuMP.add_to_expression!(objective, term) - end + # Get the set of network expansion pumps at time index `n`. + for (a, ne_pump) in ref(wm, n, :ne_pump) + # Add the cost of network expansion component `a` at time period `n`. + term = ne_pump["construction_cost"] * var(wm, n, :x_ne_pump, a) + JuMP.add_to_expression!(objective, term) end # Minimize the total cost of network expansion. @@ -180,3 +173,42 @@ function objective_owf(wm::AbstractWaterModel)::JuMP.AffExpr # Minimize the cost of network operation. return JuMP.@objective(wm.model, JuMP.MIN_SENSE, objective) end + + + +""" + objective_max_demand(wm::AbstractWaterModel)::JuMP.AffExpr +""" +function objective_max_scaled_demand(wm::AbstractWaterModel)::JuMP.AffExpr + # Get all network IDs in the multinetwork. + network_ids = sort(collect(nw_ids(wm))) + + # Find the network IDs over which the objective will be defined. + if length(network_ids) > 1 + network_ids_flow = network_ids[1:end-1] + else + network_ids_flow = network_ids + end + + # Initialize the objective expression to zero. + objective = JuMP.AffExpr(0.0) + + scale = 0.0 + for n in network_ids_flow + # Get the set of dispatchable demands at time index `n`. + dispatchable_demands = ref(wm, n, :dispatchable_demand) + + for (i, demand) in filter(x -> x.second["flow_min"] >= 0.0, dispatchable_demands) + # Add the volume delivered at demand `i` and time period `n`. + # coeff = get(demand, "priority", 1.0) * ref(wm, n, :time_step) + scale = max(scale,demand["flow_max"]) + JuMP.add_to_expression!(objective, var(wm, n, :q_demand, i)) + end + end + if(scale==0.0) + scale = 1.0 + end + println("Using this scale for water objective $(scale)*****") + # Maximize the total amount of prioritized water volume delivered. + return JuMP.@objective(wm.model, JuMP.MAX_SENSE, objective/scale) +end diff --git a/src/core/pump.jl b/src/core/pump.jl index a9659ec7..1a7e39fb 100644 --- a/src/core/pump.jl +++ b/src/core/pump.jl @@ -427,7 +427,7 @@ function _calc_pump_power_points(wm::AbstractWaterModel, nw::Int, pump_id::Int, end -function _calc_ne_pump_power_points(wm::AbstractWaterModel, nw::Int, pump_id::Int, num_points::Int) +function _calc_pump_power_points_ne(wm::AbstractWaterModel, nw::Int, pump_id::Int, num_points::Int) ne_pump = ref(wm, nw, :ne_pump, pump_id) head_curve_function = ref(wm, nw, :ne_pump, pump_id, "head_curve_function") @@ -456,16 +456,39 @@ function _calc_ne_pump_power_points(wm::AbstractWaterModel, nw::Int, pump_id::In return q_build, max.(0.0, density * gravity * inv.(eff) .* f_build) end +# pump_power_updates +function modify_pump_power_points(wm::AbstractWaterModel, nw::Int) + power_curve_unscaled = ref(wm, nw, :pump_power_curve) + bf = wm.ref[:it][wm_it_sym][:base_flow] + bm = wm.ref[:it][wm_it_sym][:base_mass] + bl = wm.ref[:it][wm_it_sym][:base_length] + bt = wm.ref[:it][wm_it_sym][:base_time] + bp = bm*(bl^2 / bt^3) -function _calc_pump_power(wm::AbstractWaterModel, nw::Int, pump_id::Int, q::Vector{Float64}) + q_points_scaled = power_curve_unscaled["q_points"] ./bf + f_points_sclaed = power_curve_unscaled["f_points"] ./bp + + return q_points_scaled, f_points_sclaed + # make the scaling consistent +end + +function _calc_pump_power_linear_coeff(wm::AbstractWaterModel, nw::Int, z::JuMP.VariableRef) + q_true, f_true = modify_pump_power_points(wm, nw) + q_array = hcat(q_true, ones(length(q_true))) + linear_coefficients = q_array \ f_true + return q -> sum(linear_coefficients .*[q, z]) +end + + +function _calc_pump_power(wm::AbstractWaterModel, nw::Int, pump_id::Int) q_true, f_true = _calc_pump_power_points(wm, nw, pump_id, 100) - return max.(Interpolations.linear_interpolation(q_true, f_true).(q), 0.0) + return q -> max.(Interpolations.linear_interpolation(q_true, f_true).(q), 0.0) end -function _calc_pump_power_ne(wm::AbstractWaterModel, nw::Int, pump_id::Int, q::Vector{Float64}) +function _calc_pump_power_ne(wm::AbstractWaterModel, nw::Int, pump_id::Int) q_true, f_true = _calc_pump_power_points_ne(wm, nw, pump_id, 100) - return max.(Interpolations.LinearInterpolation(q_true, f_true).(q), 0.0) + return q -> max.(Interpolations.linear_interpolation(q_true, f_true).(q), 0.0) end diff --git a/src/form/crd.jl b/src/form/crd.jl index 1cfa5d29..09b56ae9 100644 --- a/src/form/crd.jl +++ b/src/form/crd.jl @@ -163,59 +163,57 @@ function constraint_on_off_pump_head_gain_ne( end -function constraint_on_off_pump_power( - wm::AbstractCRDModel, - n::Int, - a::Int, - q_min_forward::Float64, -) - # Gather pump flow, power, and status variables. - q, P, z = var(wm, n, :qp_pump, a), var(wm, n, :P_pump, a), var(wm, n, :z_pump, a) - - # Compute pump flow and power partitioning. - q_lb, q_ub = q_min_forward, JuMP.upper_bound(q) - - if q_lb == 0.0 && q_ub == 0.0 - c = JuMP.@constraint(wm.model, P == 0.0) - append!(con(wm, n, :on_off_pump_power)[a], [c]) - else - f_ua = _calc_pump_power_ua(wm, n, a, [q_lb, q_ub]) - - if f_ua[1] != f_ua[2] - # Build a linear under-approximation of the power. - slope = (f_ua[2] - f_ua[1]) / (q_ub - q_lb) - power_expr = slope * (q - q_lb * z) + f_ua[1] * z - c = JuMP.@constraint(wm.model, power_expr <= P) - append!(con(wm, n, :on_off_pump_power)[a], [c]) - end - end -end - - -function constraint_on_off_pump_power_ne( - wm::AbstractCRDModel, - n::Int, - a::Int, - q_min_forward::Float64, -) - # Gather pump flow, power, and status variables. - q, P, z = var(wm, n, :qp_ne_pump, a), var(wm, n, :P_ne_pump, a), var(wm, n, :z_ne_pump, a) - - # Compute pump flow and power partitioning. - q_lb, q_ub = q_min_forward, JuMP.upper_bound(q) - - if q_lb == 0.0 && q_ub == 0.0 - c = JuMP.@constraint(wm.model, P == 0.0) - append!(con(wm, n, :on_off_pump_power_ne)[a], [c]) - else - f_ua = _calc_pump_power_ua_ne(wm, n, a, [q_lb, q_ub]) - - if f_ua[1] != f_ua[2] - # Build a linear under-approximation of the power. - slope = (f_ua[2] - f_ua[1]) / (q_ub - q_lb) - power_expr = slope * (q - q_lb * z) + f_ua[1] * z - c = JuMP.@constraint(wm.model, power_expr <= P) - append!(con(wm, n, :on_off_pump_power_ne)[a], [c]) - end - end -end +# function constraint_on_off_pump_power( +# wm::AbstractCRDModel, +# n::Int, +# a::Int, +# q_min_forward::Float64, +# ) +# # Gather pump flow, power, and status variables. +# q, P, z = var(wm, n, :qp_pump, a), var(wm, n, :P_pump, a), var(wm, n, :z_pump, a) +# +# # Compute pump flow and power partitioning. +# q_lb, q_ub = q_min_forward, JuMP.upper_bound(q) +# +# if q_lb == 0.0 && q_ub == 0.0 +# c = JuMP.@constraint(wm.model, P == 0.0) +# append!(con(wm, n, :on_off_pump_power)[a], [c]) +# else +# f_ua = _calc_pump_power_ua(wm, n, a, [q_lb, q_ub]) +# if f_ua[1] != f_ua[2] +# # Build a linear under-approximation of the power. +# slope = (f_ua[2] - f_ua[1]) / (q_ub - q_lb) +# power_expr = slope * (q - q_lb * z) + f_ua[1] * z +# c = JuMP.@constraint(wm.model, power_expr <= P) +# append!(con(wm, n, :on_off_pump_power)[a], [c]) +# end +# end +# end +# +# +# function constraint_on_off_pump_power_ne( +# wm::AbstractCRDModel, +# n::Int, +# a::Int, +# q_min_forward::Float64, +# ) +# # Gather pump flow, power, and status variables. +# q, P, z = var(wm, n, :qp_ne_pump, a), var(wm, n, :P_ne_pump, a), var(wm, n, :z_ne_pump, a) +# +# # Compute pump flow and power partitioning. +# q_lb, q_ub = q_min_forward, JuMP.upper_bound(q) +# +# if q_lb == 0.0 && q_ub == 0.0 +# c = JuMP.@constraint(wm.model, P == 0.0) +# append!(con(wm, n, :on_off_pump_power_ne)[a], [c]) +# else +# f_ua = _calc_pump_power_ua_ne(wm, n, a, [q_lb, q_ub]) +# if f_ua[1] != f_ua[2] +# # Build a linear under-approximation of the power. +# slope = (f_ua[2] - f_ua[1]) / (q_ub - q_lb) +# power_expr = slope * (q - q_lb * z) + f_ua[1] * z +# c = JuMP.@constraint(wm.model, power_expr <= P) +# append!(con(wm, n, :on_off_pump_power_ne)[a], [c]) +# end +# end +# end diff --git a/src/form/nc.jl b/src/form/nc.jl index 5a3b687f..20a5bc71 100644 --- a/src/form/nc.jl +++ b/src/form/nc.jl @@ -190,8 +190,12 @@ function constraint_pipe_head_loss( q, h_i, h_j = var(wm, n, :q_pipe, a), var(wm, n, :h, node_fr), var(wm, n, :h, node_to) # Add nonconvex constraint for the head loss relationship. - c_1 = JuMP.@NLconstraint(wm.model, r * head_loss(q) <= (h_i - h_j) / L) - c_2 = JuMP.@NLconstraint(wm.model, r * head_loss(q) >= (h_i - h_j) / L) + head_loss_form = wm.ref[:it][wm_it_sym][:head_loss] + p = uppercase(head_loss_form) == "H-W" ? 1.852 : 2.0 + alpha = (p-1) + h_l = JuMP.@expression(wm.model, q*(abs(q)^alpha)) + c_1 = JuMP.@NLconstraint(wm.model, r * h_l <= (h_i - h_j) / L) + c_2 = JuMP.@NLconstraint(wm.model, r * h_l >= (h_i - h_j) / L) # Append the :pipe_head_loss constraint array. append!(con(wm, n, :pipe_head_loss)[a], [c_1, c_2]) @@ -610,16 +614,19 @@ function constraint_on_off_pump_power( a::Int, q_min_forward::Float64, ) - # Gather pump flow, scaled power, and status variables. - q, P, z = var(wm, n, :q_pump, a), var(wm, n, :P_pump, a), var(wm, n, :z_pump, a) - - # Add constraint equating power with respect to the power curve. - power_qa = _calc_pump_power_quadratic_approximation(wm, n, a, z) - c_1 = JuMP.@constraint(wm.model, power_qa(q) <= P) - c_2 = JuMP.@constraint(wm.model, power_qa(q) >= P) - - # Append the :on_off_pump_power constraint array. - append!(con(wm, n, :on_off_pump_power)[a], [c_1, c_2]) + # Gather pump flow, power, and status variables. + q = var(wm, n, :q_pump, a) + P = var(wm, n, :P_pump, a) + z = var(wm, n, :z_pump, a) + + # Add constraint equating power with respect to the power curve. + power_la = _calc_pump_power_linear_coeff(wm, n, z) + # power_la = _calc_pump_power_quadratic_approximation(wm, n, a, z) + c_1 = JuMP.@constraint(wm.model, power_la(q) <= P) + c_2 = JuMP.@constraint(wm.model, power_la(q) >= P) + + # Append the :on_off_pump_power constraint array. + append!(con(wm, n, :on_off_pump_power)[a], [c_1, c_2]) end @@ -629,16 +636,19 @@ function constraint_on_off_pump_power_ne( a::Int, q_min_forward::Float64, ) - # Gather pump flow, scaled power, and status variables. - q, P, z = var(wm, n, :q_ne_pump, a), var(wm, n, :P_ne_pump, a), var(wm, n, :z_ne_pump, a) - - # Add constraint equating power with respect to the power curve. - power_qa = _calc_pump_power_quadratic_approximation_ne(wm, n, a, z) - c_1 = JuMP.@constraint(wm.model, power_qa(q) <= P) - c_2 = JuMP.@constraint(wm.model, power_qa(q) >= P) - - # Append the :on_off_pump_power constraint array. - append!(con(wm, n, :on_off_pump_power_ne)[a], [c_1, c_2]) + # Gather pump flow, power, and status variables. + q = var(wm, n, :q_ne_pump, a) + P = var(wm, n, :P_ne_pump, a) + z = var(wm, n, :z_ne_pump, a) + + # Add constraint equating power with respect to the power curve. + # power_la = _calc_pump_power_quadratic_approximation_ne(wm, n, a, z) + power_la = _calc_pump_power_linear_coeff(wm, n, z) + c_1 = JuMP.@constraint(wm.model, power_la(q) <= P) + c_2 = JuMP.@constraint(wm.model, power_la(q) >= P) + + # Append the :on_off_pump_power constraint array. + append!(con(wm, n, :on_off_pump_power_ne)[a], [c_1, c_2]) end @@ -848,6 +858,7 @@ function constraint_short_pipe_flow_ne( q_max_reverse::Float64, q_min_forward::Float64, ) +# println("Using NC ******* short pipe *********") # Get flow and status variables for the short pipe. q = var(wm, n, :q_ne_short_pipe, a) z = var(wm, n, :z_ne_short_pipe, a) diff --git a/src/form/ncd.jl b/src/form/ncd.jl index 648286d6..f6881a29 100644 --- a/src/form/ncd.jl +++ b/src/form/ncd.jl @@ -343,13 +343,15 @@ function constraint_pipe_head_loss( ) # Add constraints for positive flow and head difference. qp, dhp = var(wm, n, :qp_pipe, a), var(wm, n, :dhp_pipe, a) - c_1 = JuMP.@NLconstraint(wm.model, r * head_loss(qp) <= dhp / L) - c_2 = JuMP.@NLconstraint(wm.model, r * head_loss(qp) >= dhp / L) + head_loss_form = wm.ref[:it][wm_it_sym][:head_loss] + p = uppercase(head_loss_form) == "H-W" ? 1.852 : 2.0 + c_1 = JuMP.@NLconstraint(wm.model, r * (qp^p) <= dhp / L) + c_2 = JuMP.@NLconstraint(wm.model, r * (qp^p) >= dhp / L) # Add constraints for negative flow and head difference. qn, dhn = var(wm, n, :qn_pipe, a), var(wm, n, :dhn_pipe, a) - c_3 = JuMP.@NLconstraint(wm.model, r * head_loss(qn) <= dhn / L) - c_4 = JuMP.@NLconstraint(wm.model, r * head_loss(qn) >= dhn / L) + c_3 = JuMP.@NLconstraint(wm.model, r * (qn^p) <= dhn / L) + c_4 = JuMP.@NLconstraint(wm.model, r * (qn^p) >= dhn / L) # Append the :pipe_head_loss constraint array. append!(con(wm, n, :pipe_head_loss)[a], [c_1, c_2, c_3, c_4]) @@ -735,9 +737,10 @@ function constraint_on_off_pump_power( z = var(wm, n, :z_pump, a) # Add constraint equating power with respect to the power curve. - power_qa = _calc_pump_power_quadratic_approximation(wm, n, a, z) - c_1 = JuMP.@constraint(wm.model, power_qa(q) <= P) - c_2 = JuMP.@constraint(wm.model, power_qa(q) >= P) + power_la = _calc_pump_power_linear_coeff(wm, n, z) + # power_la = _calc_pump_power_quadratic_approximation(wm, n, a, z) + c_1 = JuMP.@constraint(wm.model, power_la(q) <= P) + c_2 = JuMP.@constraint(wm.model, power_la(q) >= P) # Append the :on_off_pump_power constraint array. append!(con(wm, n, :on_off_pump_power)[a], [c_1, c_2]) @@ -756,9 +759,10 @@ function constraint_on_off_pump_power_ne( z = var(wm, n, :z_ne_pump, a) # Add constraint equating power with respect to the power curve. - power_qa = _calc_pump_power_quadratic_approximation_ne(wm, n, a, z) - c_1 = JuMP.@constraint(wm.model, power_qa(q) <= P) - c_2 = JuMP.@constraint(wm.model, power_qa(q) >= P) + # power_la = _calc_pump_power_quadratic_approximation_ne(wm, n, a, z) + power_la = _calc_pump_power_linear_coeff(wm, n, z) + c_1 = JuMP.@constraint(wm.model, power_la(q) <= P) + c_2 = JuMP.@constraint(wm.model, power_la(q) >= P) # Append the :on_off_pump_power constraint array. append!(con(wm, n, :on_off_pump_power_ne)[a], [c_1, c_2]) @@ -899,6 +903,7 @@ function constraint_short_pipe_flow_ne( q_max_reverse::Float64, q_min_forward::Float64, ) +# println("Using NCD ******* short pipe *********") # Get expansion short pipe flow, direction, and status variables. qp, qn = var(wm, n, :qp_ne_short_pipe, a), var(wm, n, :qn_ne_short_pipe, a) y, z = var(wm, n, :y_ne_short_pipe, a), var(wm, n, :z_ne_short_pipe, a) diff --git a/src/form/pwlrd.jl b/src/form/pwlrd.jl index a00ca6d2..a3cdeadc 100644 --- a/src/form/pwlrd.jl +++ b/src/form/pwlrd.jl @@ -43,7 +43,7 @@ function variable_flow(wm::AbstractPWLRDModel; nw::Int=nw_id_default, bounded::B # Create directed head difference (`dhp` and `dhn`) variables for each component. _variable_component_head_difference(wm, name; nw = nw, bounded = bounded, report = report) end - + # Create variables involved in convex combination constraints for pumps. pump_partition = Dict{Int, Vector{Float64}}(a => pump["flow_partition"] for (a, pump) in ref(wm, nw, :pump)) @@ -58,6 +58,20 @@ function variable_flow(wm::AbstractPWLRDModel; nw::Int=nw_id_default, bounded::B base_name = "$(nw)_x", binary = true, start = comp_start_value(ref(wm, nw, :pump, a), "x_start")) + # Create variables involved in convex combination constraints for expansion pumps. + pump_partition = Dict{Int, Vector{Float64}}(a => + pump["flow_partition"] for (a, pump) in ref(wm, nw, :ne_pump)) + + var(wm, nw)[:lambda_ne_pump] = JuMP.@variable(wm.model, + [a in ids(wm, nw, :ne_pump), k in 1:length(pump_partition[a])], + base_name = "$(nw)_ne_lambda", lower_bound = 0.0, upper_bound = 1.0, + start = comp_start_value(ref(wm, nw, :ne_pump, a), "lambda_start", k)) + + var(wm, nw)[:x_con_ne_pump] = JuMP.@variable(wm.model, + [a in ids(wm, nw, :ne_pump), k in 1:length(pump_partition[a]) - 1], + base_name = "$(nw)_ne_x_con", binary = true, + start = comp_start_value(ref(wm, nw, :ne_pump, a), "x_con_ne_start")) + # Create variables involved in convex combination constraints for pipes. pipe_partition_p = Dict{Int, Vector{Float64}}(a => get_pipe_flow_partition_positive(pipe) @@ -110,7 +124,7 @@ function variable_flow(wm::AbstractPWLRDModel; nw::Int=nw_id_default, bounded::B [a in ids(wm, nw, :des_pipe), k in 1:length(des_pipe_partition_n[a])], base_name = "$(nw)_lambda_n", lower_bound = 0.0, upper_bound = 1.0, start = comp_start_value(ref(wm, nw, :des_pipe, a), "lambda_n_start", k)) - + var(wm, nw)[:x_n_des_pipe] = JuMP.@variable(wm.model, [a in ids(wm, nw, :des_pipe), k in 1:length(des_pipe_partition_n[a])-1], base_name = "$(nw)_x_n", binary = true, @@ -205,7 +219,7 @@ function constraint_pipe_head_loss( for flow_value in filter(x -> x > 0.0, partition_p) # Add head loss outer (i.e., lower) approximations. lhs_p = r * _calc_head_loss_oa(qp, y, flow_value, exponent) - + if minimum(abs.(lhs_p.terms.vals)) >= 1.0e-4 scalar = _get_scaling_factor(vcat(lhs_p.terms.vals, [1.0 / L])) c_7 = JuMP.@constraint(wm.model, scalar * lhs_p <= scalar * dhp / L) @@ -243,7 +257,7 @@ function constraint_pipe_head_loss( c_13 = JuMP.@constraint(wm.model, lambda_n[a, bn_range[end]] <= x_n[a, bn_range_m1[end]]) append!(con(wm, n, :pipe_head_loss, a), [c_11, c_12, c_13]) end - + # Add a constraint that upper-bounds the head loss variable. if maximum(partition_n) != 0.0 f_n = r .* partition_n.^exponent @@ -412,7 +426,7 @@ function constraint_on_off_des_pipe_head_loss( lhs_n_2 = r * _calc_head_loss_oa(qn, z, flow_value, exponent) scalar = _get_scaling_factor(vcat(lhs_n_2.terms.vals, [1.0 / L])) c_17 = JuMP.@constraint(wm.model, scalar * lhs_n_2 <= scalar * dhn / L) - + append!(con(wm, n, :on_off_des_pipe_head_loss, a), [c_16, c_17]) end @@ -428,7 +442,7 @@ function constraint_on_off_des_pipe_head_loss( c_19 = JuMP.@constraint(wm.model, x_p_sum + x_n_sum == z) lambda_p_sum = sum(lambda_p[a, k] for k in bp_range) - lambda_n_sum = sum(lambda_n[a, k] for k in bn_range) + lambda_n_sum = sum(lambda_n[a, k] for k in bn_range) c_20 = JuMP.@constraint(wm.model, lambda_p_sum + lambda_n_sum == z) append!(con(wm, n, :on_off_des_pipe_head_loss, a), [c_19, c_20]) end @@ -523,50 +537,135 @@ end """ - constraint_on_off_pump_power( + constraint_on_off_pump_head_gain( wm::AbstractPWLRDModel, n::Int, a::Int, + node_fr::Int, + node_to::Int, q_min_forward::Float64 ) -Adds constraints that model the pump's power consumption, if operating, as -being bounded between piecewise-linear inner and outer approximations of a -(potentially) nonlinear function of a more accurate model. If the pump is -inactive, the power is restricted to a value of zero. Here, `wm` is the -WaterModels object, `n` is the subnetwork (or time) index that is considered, -`a` is the index of the pump, and `q_min_forward` is the _minimum_ (positive) -amount of flow when the pump is active. +Adds constraints that model the pump's head gain, if operating, via outer +and piecewise-linear inner approximations of the nonlinear function of flow +rate. If the pump is inactive, the head gain is restricted to a value of zero. +Here, `wm` is the WaterModels object, `n` is the subnetwork (or time) index +that is considered, `a` is the index of the pump, `node_fr` is the index of the +tail node of the pump, `node_to` is the index of the head node of the pump, and +`q_min_forward` is the _minimum_ (positive) amount of flow when the pump is +active. Head gain is assumed to be nonnegative and is directed from `node_fr` +to `node_to`. """ -function constraint_on_off_pump_power( +function constraint_on_off_pump_head_gain_ne( wm::AbstractPWLRDModel, n::Int, a::Int, + node_fr::Int, + node_to::Int, q_min_forward::Float64 ) - # Gather pump power and status variables. - P, lambda = var(wm, n, :P_pump, a), var(wm, n, :lambda_pump) + # Gather pump flow, head gain, and status variables. + qp, g, z = var(wm, n, :qp_ne_pump, a), var(wm, n, :g_ne_pump, a), var(wm, n, :z_ne_pump, a) + + # Gather convex combination variables. + lambda, x = var(wm, n, :lambda_ne_pump), var(wm, n, :x_con_ne_pump) # Get the corresponding flow partitioning. - @assert haskey(ref(wm, n, :pump, a), "flow_partition") - partition = ref(wm, n, :pump, a, "flow_partition") - bp_range = 1:length(partition) - - # Add a constraint that lower-bounds the power variable. - if minimum(partition) == 0.0 && maximum(partition) == 0.0 - c_1 = JuMP.@constraint(wm.model, P == 0.0) - append!(con(wm, n, :on_off_pump_power)[a], [c_1]) - else - f_ua = _calc_pump_power_ua(wm, n, a, partition) - power_lb_expr = sum(f_ua[k] * lambda[a, k] for k in bp_range) - c_1 = JuMP.@constraint(wm.model, power_lb_expr <= P) + @assert haskey(ref(wm, n, :ne_pump, a), "flow_partition") + partition = ref(wm, n, :ne_pump, a, "flow_partition") + bp_range, bp_range_m1 = 1:length(partition), 1:length(partition)-1 + + # Add the required SOS constraints. + c_1 = JuMP.@constraint(wm.model, sum(lambda[a, k] for k in bp_range) == z) + c_2 = JuMP.@constraint(wm.model, sum(x[a, k] for k in bp_range_m1) == z) + c_3 = JuMP.@constraint(wm.model, lambda[a, 1] <= x[a, 1]) + c_4 = JuMP.@constraint(wm.model, lambda[a, bp_range[end]] <= x[a, bp_range_m1[end]]) + + # Add a constraint for the flow piecewise approximation. + qp_lhs = sum(partition[k] * lambda[a, k] for k in bp_range) + scalar = _get_scaling_factor(vcat(qp_lhs.terms.vals, [1.0])) + c_5 = JuMP.@constraint(wm.model, scalar * qp_lhs == scalar * qp) + + # Get pump head curve function and its derivative. + head_curve_function = ref(wm, n, :ne_pump, a, "head_curve_function") + head_curve_derivative = ref(wm, n, :ne_pump, a, "head_curve_derivative") + + # Add a constraint that lower-bounds the head gain variable. + f_all = map(x -> x = x < 0.0 ? 0.0 : x, head_curve_function.(partition)) + gain_lb_expr = sum(f_all[k] .* lambda[a, k] for k in bp_range) + scalar = _get_scaling_factor(vcat(gain_lb_expr.terms.vals, [1.0])) + c_6 = JuMP.@constraint(wm.model, scalar * gain_lb_expr <= scalar * g) - # Add a constraint that upper-bounds the power variable. - f_oa = _calc_pump_power_oa(wm, n, a, partition) - power_ub_expr = sum(f_oa[k] * lambda[a, k] for k in bp_range) - c_2 = JuMP.@constraint(wm.model, P <= power_ub_expr) + # Append the constraint array. + append!(con(wm, n, :on_off_pump_head_gain_ne, a), [c_1, c_2, c_3, c_4, c_5, c_6]) + + for flow_value in filter(q_val -> q_val > 0.0, partition) + # Add head gain outer (i.e., upper) approximations. + f = head_curve_function(flow_value) + df = head_curve_derivative(flow_value) + + if abs(df) >= 1.0e-4 # Only add an outer-approximation if the derivative isn't too small. + f_rhs = f * z + df * (qp - flow_value * z) + scalar = _get_scaling_factor(vcat(f_rhs.terms.vals, [1.0])) + c = JuMP.@constraint(wm.model, scalar * g <= scalar * f_rhs) + append!(con(wm, n, :on_off_pump_head_gain_ne, a), [c]) + end + end - # Append the :on_off_pump_power constraint array. - append!(con(wm, n, :on_off_pump_power)[a], [c_1, c_2]) + for k in 2:length(partition)-1 + # Add the adjacency constraints for piecewise variables. + adjacency = x[a, k-1] + x[a, k] + c = JuMP.@constraint(wm.model, lambda[a, k] <= adjacency) + append!(con(wm, n, :on_off_pump_head_gain_ne, a), [c]) end end + + +""" + constraint_on_off_pump_power( + wm::AbstractPWLRDModel, + n::Int, + a::Int, + q_min_forward::Float64 + ) + +Adds constraints that model the pump's power consumption, if operating, as +being bounded between piecewise-linear inner and outer approximations of a +(potentially) nonlinear function of a more accurate model. If the pump is +inactive, the power is restricted to a value of zero. Here, `wm` is the +WaterModels object, `n` is the subnetwork (or time) index that is considered, +`a` is the index of the pump, and `q_min_forward` is the _minimum_ (positive) +amount of flow when the pump is active. +""" +# function constraint_on_off_pump_power( +# wm::AbstractPWLRDModel, +# n::Int, +# a::Int, +# q_min_forward::Float64 +# ) +# # Gather pump power and status variables. +# P, lambda = var(wm, n, :P_pump, a), var(wm, n, :lambda_pump) +# +# # Get the corresponding flow partitioning. +# @assert haskey(ref(wm, n, :pump, a), "flow_partition") +# partition = ref(wm, n, :pump, a, "flow_partition") +# bp_range = 1:length(partition) +# +# # Add a constraint that lower-bounds the power variable. +# if minimum(partition) == 0.0 && maximum(partition) == 0.0 +# c_1 = JuMP.@constraint(wm.model, P == 0.0) +# append!(con(wm, n, :on_off_pump_power)[a], [c_1]) +# else +# f_ua = _calc_pump_power_ua(wm, n, a, partition) +# power_lb_expr = sum(f_ua[k] * lambda[a, k] for k in bp_range) +# c_1 = JuMP.@constraint(wm.model, power_lb_expr <= P) +# +# # Add a constraint that upper-bounds the power variable. +# f_oa = _calc_pump_power_oa(wm, n, a, partition) +# power_ub_expr = sum(f_oa[k] * lambda[a, k] for k in bp_range) +# c_2 = JuMP.@constraint(wm.model, P <= power_ub_expr) +# +# # Append the :on_off_pump_power constraint array. +# append!(con(wm, n, :on_off_pump_power)[a], [c_1, c_2]) +# end +# end diff --git a/src/form/scip_constraints_variables.jl b/src/form/scip_constraints_variables.jl new file mode 100644 index 00000000..6178f038 --- /dev/null +++ b/src/form/scip_constraints_variables.jl @@ -0,0 +1,358 @@ +# This file is added as a temporary fix to work with scip. +# In future, the pipe head_loss function needs to be changed + +function constraint_pipe_head_loss_scip_lp( + wm::AbstractWaterModel, + a::Int; + nw::Int = nw_id_default, + kwargs..., +) + node_fr = ref(wm, nw, :pipe, a, "node_fr") + node_to = ref(wm, nw, :pipe, a, "node_to") + exponent = _get_exponent_from_head_loss_form(wm.ref[:it][wm_it_sym][:head_loss]) + L = ref(wm, nw, :pipe, a, "length") + + wm_data = get_wm_data(wm.data) + head_loss, viscosity = wm_data["head_loss"], wm_data["viscosity"] + base_length = wm_data["per_unit"] ? wm_data["base_length"] : 1.0 + base_mass = wm_data["per_unit"] ? wm_data["base_mass"] : 1.0 + base_time = wm_data["per_unit"] ? wm_data["base_time"] : 1.0 + + r = _calc_pipe_resistance( + ref(wm, nw, :pipe, a), + head_loss, + viscosity, + base_length, + base_mass, + base_time, + ) + q_max_reverse = min(get(ref(wm, nw, :pipe, a), "flow_max_reverse", 0.0), 0.0) + q_min_forward = max(get(ref(wm, nw, :pipe, a), "flow_min_forward", 0.0), 0.0) + + _initialize_con_dict(wm, :pipe_head_loss, nw = nw, is_array = true) + con(wm, nw, :pipe_head_loss)[a] = Array{JuMP.ConstraintRef}([]) + constraint_pipe_head_loss_scip_lp( + wm, + nw, + a, + node_fr, + node_to, + exponent, + L, + r, + q_max_reverse, + q_min_forward, + ) +end + + + +function constraint_pipe_head_loss_scip_lp( + wm::AbstractNCDModel, + n::Int, + a::Int, + node_fr::Int, + node_to::Int, + exponent::Float64, + L::Float64, + r::Float64, + q_max_reverse::Float64, + q_min_forward::Float64, +) +# Get the variable for flow directionality. +y = var(wm, n, :y_pipe, a) + +# Get variables for positive flow and head difference. +qp, dhp = var(wm, n, :qp_pipe, a), var(wm, n, :dhp_pipe, a) + +# Get positively-directed convex combination variables. +lambda_p, x_p = var(wm, n, :lambda_p_pipe), var(wm, n, :x_p_pipe) + +# Get the corresponding positive flow partitioning. +partition_p = get_pipe_flow_partition_positive(ref(wm, n, :pipe, a)) +bp_range, bp_range_m1 = 1:length(partition_p), 1:length(partition_p)-1 + +# Add constraints for the positive flow piecewise approximation. +c_1 = JuMP.@constraint(wm.model, sum(lambda_p[a, k] for k in bp_range) == y) +qp_sum = sum(partition_p[k] * lambda_p[a, k] for k in bp_range) +scalar = _get_scaling_factor(vcat(qp_sum.terms.vals, [1.0])) +c_2 = JuMP.@constraint(wm.model, scalar * qp_sum == scalar * qp) +append!(con(wm, n, :pipe_head_loss, a), [c_1, c_2]) + +if length(partition_p) > 1 + # If there are multiple points, constrain the convex combination. + c_3 = JuMP.@constraint(wm.model, sum(x_p[a, k] for k in bp_range_m1) == y) + c_4 = JuMP.@constraint(wm.model, lambda_p[a, 1] <= x_p[a, 1]) + c_5 = JuMP.@constraint(wm.model, lambda_p[a, bp_range[end]] <= x_p[a, bp_range_m1[end]]) + append!(con(wm, n, :pipe_head_loss, a), [c_3, c_4, c_5]) +end + +# Add a constraint that upper-bounds the head loss variable. +if maximum(partition_p) != 0.0 + f_p = r * partition_p.^exponent + f_p_ub_expr = sum(f_p[k] * lambda_p[a, k] for k in bp_range) + scalar = _get_scaling_factor(vcat(f_p_ub_expr.terms.vals, [1.0 / L])) + c_6 = JuMP.@constraint(wm.model, scalar * dhp / L <= scalar * f_p_ub_expr) + append!(con(wm, n, :pipe_head_loss, a), [c_6]) +else + c_6 = JuMP.@constraint(wm.model, dhp == 0.0) + append!(con(wm, n, :pipe_head_loss, a), [c_6]) +end + +# Loop over consequential points (i.e., those that have nonzero head loss). +for flow_value in filter(x -> x > 0.0, partition_p) + # Add head loss outer (i.e., lower) approximations. + lhs_p = r * _calc_head_loss_oa(qp, y, flow_value, exponent) + + if minimum(abs.(lhs_p.terms.vals)) >= 1.0e-4 + scalar = _get_scaling_factor(vcat(lhs_p.terms.vals, [1.0 / L])) + c_7 = JuMP.@constraint(wm.model, scalar * lhs_p <= scalar * dhp / L) + append!(con(wm, n, :pipe_head_loss, a), [c_7]) + end +end + +for k in 2:length(partition_p)-1 + # Add the adjacency constraints for piecewise variables. + c_8 = JuMP.@constraint(wm.model, lambda_p[a, k] <= x_p[a, k-1] + x_p[a, k]) + append!(con(wm, n, :pipe_head_loss, a), [c_8]) +end + +# Get variables for negative flow and head difference. +qn, dhn = var(wm, n, :qn_pipe, a), var(wm, n, :dhn_pipe, a) + +# Get negatively-directed convex combination variable. +lambda_n, x_n = var(wm, n, :lambda_n_pipe), var(wm, n, :x_n_pipe) + +# Get the corresponding negative flow partitioning (negated). +partition_n = sort(-get_pipe_flow_partition_negative(ref(wm, n, :pipe, a))) +bn_range, bn_range_m1 = 1:length(partition_n), 1:length(partition_n)-1 + +# Add constraints for the negative flow piecewise approximation. +c_9 = JuMP.@constraint(wm.model, sum(lambda_n[a, k] for k in bn_range) == 1.0 - y) +qn_sum = sum(partition_n[k] * lambda_n[a, k] for k in bn_range) +scalar = _get_scaling_factor(vcat(qn_sum.terms.vals, [1.0])) +c_10 = JuMP.@constraint(wm.model, scalar * qn_sum == scalar * qn) +append!(con(wm, n, :pipe_head_loss, a), [c_9, c_10]) + +if length(partition_n) > 1 + # If there are multiple points, constrain the convex combination. + c_11 = JuMP.@constraint(wm.model, sum(x_n[a, k] for k in bn_range_m1) == 1.0 - y) + c_12 = JuMP.@constraint(wm.model, lambda_n[a, 1] <= x_n[a, 1]) + c_13 = JuMP.@constraint(wm.model, lambda_n[a, bn_range[end]] <= x_n[a, bn_range_m1[end]]) + append!(con(wm, n, :pipe_head_loss, a), [c_11, c_12, c_13]) +end + +# Add a constraint that upper-bounds the head loss variable. +if maximum(partition_n) != 0.0 + f_n = r .* partition_n.^exponent + f_n_ub_expr = sum(f_n[k] * lambda_n[a, k] for k in bn_range) + scalar = _get_scaling_factor(vcat(f_n_ub_expr.terms.vals, [1.0 / L])) + c_14 = JuMP.@constraint(wm.model, scalar * dhn / L <= scalar * f_n_ub_expr) + append!(con(wm, n, :pipe_head_loss, a), [c_14]) +else + c_14 = JuMP.@constraint(wm.model, dhn == 0.0) + append!(con(wm, n, :pipe_head_loss, a), [c_14]) +end + +# Loop over consequential points (i.e., those that have nonzero head loss). +for flow_value in filter(x -> x > 0.0, partition_n) + # Add head loss outer (i.e., lower) approximations. + lhs_n = r * _calc_head_loss_oa(qn, 1.0 - y, flow_value, exponent) + + if minimum(abs.(lhs_n.terms.vals)) >= 1.0e-4 + scalar = _get_scaling_factor(vcat(lhs_n.terms.vals, [1.0 / L])) + c_15 = JuMP.@constraint(wm.model, scalar * lhs_n <= scalar * dhn / L) + append!(con(wm, n, :pipe_head_loss, a), [c_15]) + end +end + +for k in 2:length(partition_n)-1 + # Add the adjacency constraints for piecewise variables. + c_16 = JuMP.@constraint(wm.model, lambda_n[a, k] <= x_n[a, k-1] + x_n[a, k]) + append!(con(wm, n, :pipe_head_loss, a), [c_16]) +end +end + + + +# function variable_flow_scip( +# wm::AbstractNCDModel; +# nw::Int = nw_id_default, +# bounded::Bool = true, +# report::Bool = true, +# ) +# for name in _LINK_COMPONENTS +# # Create directed flow (`qp` and `qn`) variables for each component. +# _variable_component_flow(wm, name; nw = nw, bounded = bounded, report = report) +# +# # Create directed flow binary direction variables (`y`) for each component. +# _variable_component_direction(wm, name; nw = nw, report = report) +# end +# +# +# # Create variables involved in convex combination constraints for pipes. +# pipe_partition_p = Dict{Int, Vector{Float64}}(a => +# get_pipe_flow_partition_positive(pipe) +# for (a, pipe) in ref(wm, nw, :pipe)) +# +# var(wm, nw)[:lambda_p_pipe] = JuMP.@variable(wm.model, +# [a in ids(wm, nw, :pipe), k in 1:length(pipe_partition_p[a])], +# base_name = "$(nw)_lambda_p", lower_bound = 0.0, upper_bound = 1.0, +# start = comp_start_value(ref(wm, nw, :pipe, a), "lambda_p_start", k)) +# +# var(wm, nw)[:x_p_pipe] = JuMP.@variable(wm.model, +# [a in ids(wm, nw, :pipe), k in 1:length(pipe_partition_p[a])-1], +# base_name = "$(nw)_x_p", binary = true, +# start = comp_start_value(ref(wm, nw, :pipe, a), "x_p_start")) +# +# pipe_partition_n = Dict{Int, Vector{Float64}}(a => +# sort(-get_pipe_flow_partition_negative(pipe)) +# for (a, pipe) in ref(wm, nw, :pipe)) +# +# var(wm, nw)[:lambda_n_pipe] = JuMP.@variable(wm.model, +# [a in ids(wm, nw, :pipe), k in 1:length(pipe_partition_n[a])], +# base_name = "$(nw)_lambda_n", lower_bound = 0.0, upper_bound = 1.0, +# start = comp_start_value(ref(wm, nw, :pipe, a), "lambda_n_start", k)) +# +# var(wm, nw)[:x_n_pipe] = JuMP.@variable(wm.model, +# [a in ids(wm, nw, :pipe), k in 1:length(pipe_partition_n[a])-1], +# base_name = "$(nw)_x_n", binary = true, +# start = comp_start_value(ref(wm, nw, :pipe, a), "x_n_start")) +# +# for name in ["des_pipe", "pipe"] +# # Create directed head difference (`dhp` and `dhn`) variables for each component. +# _variable_component_head_difference( +# wm, +# name; +# nw = nw, +# bounded = bounded, +# report = report, +# ) +# end +# end + +# function constraint_pipe_head_loss_scip( +# wm::AbstractNCDModel, +# n::Int, +# a::Int, +# node_fr::Int, +# node_to::Int, +# exponent::Float64, +# L::Float64, +# r::Float64, +# q_max_reverse::Float64, +# q_min_forward::Float64 +# ) +# # Get the variable for flow directionality. +# y = var(wm, n, :y_pipe, a) +# +# # Get variables for positive flow and head difference. +# qp, dhp = var(wm, n, :qp_pipe, a), var(wm, n, :dhp_pipe, a) +# +# # Get positively-directed convex combination variables. +# lambda_p, x_p = var(wm, n, :lambda_p_pipe), var(wm, n, :x_p_pipe) +# +# # Get the corresponding positive flow partitioning. +# partition_p = get_pipe_flow_partition_positive(ref(wm, n, :pipe, a)) +# bp_range, bp_range_m1 = 1:length(partition_p), 1:length(partition_p)-1 +# +# # Add constraints for the positive flow piecewise approximation. +# c_1 = JuMP.@constraint(wm.model, sum(lambda_p[a, k] for k in bp_range) == y) +# qp_sum = sum(partition_p[k] * lambda_p[a, k] for k in bp_range) +# scalar = _get_scaling_factor(vcat(qp_sum.terms.vals, [1.0])) +# c_2 = JuMP.@constraint(wm.model, scalar * qp_sum == scalar * qp) +# append!(con(wm, n, :pipe_head_loss, a), [c_1, c_2]) +# +# if length(partition_p) > 1 +# # If there are multiple points, constrain the convex combination. +# c_3 = JuMP.@constraint(wm.model, sum(x_p[a, k] for k in bp_range_m1) == y) +# c_4 = JuMP.@constraint(wm.model, lambda_p[a, 1] <= x_p[a, 1]) +# c_5 = JuMP.@constraint(wm.model, lambda_p[a, bp_range[end]] <= x_p[a, bp_range_m1[end]]) +# append!(con(wm, n, :pipe_head_loss, a), [c_3, c_4, c_5]) +# end +# +# # Add a constraint that upper-bounds the head loss variable. +# if maximum(partition_p) != 0.0 +# f_p = r * partition_p.^exponent +# f_p_ub_expr = sum(f_p[k] * lambda_p[a, k] for k in bp_range) +# scalar = _get_scaling_factor(vcat(f_p_ub_expr.terms.vals, [1.0 / L])) +# c_6 = JuMP.@constraint(wm.model, scalar * dhp / L <= scalar * f_p_ub_expr) +# append!(con(wm, n, :pipe_head_loss, a), [c_6]) +# else +# c_6 = JuMP.@constraint(wm.model, dhp == 0.0) +# append!(con(wm, n, :pipe_head_loss, a), [c_6]) +# end +# +# # Loop over consequential points (i.e., those that have nonzero head loss). +# for flow_value in filter(x -> x > 0.0, partition_p) +# # Add head loss outer (i.e., lower) approximations. +# lhs_p = r * _calc_head_loss_oa(qp, y, flow_value, exponent) +# +# if minimum(abs.(lhs_p.terms.vals)) >= 1.0e-4 +# scalar = _get_scaling_factor(vcat(lhs_p.terms.vals, [1.0 / L])) +# c_7 = JuMP.@constraint(wm.model, scalar * lhs_p <= scalar * dhp / L) +# append!(con(wm, n, :pipe_head_loss, a), [c_7]) +# end +# end +# +# for k in 2:length(partition_p)-1 +# # Add the adjacency constraints for piecewise variables. +# c_8 = JuMP.@constraint(wm.model, lambda_p[a, k] <= x_p[a, k-1] + x_p[a, k]) +# append!(con(wm, n, :pipe_head_loss, a), [c_8]) +# end +# +# # Get variables for negative flow and head difference. +# qn, dhn = var(wm, n, :qn_pipe, a), var(wm, n, :dhn_pipe, a) +# +# # Get negatively-directed convex combination variable. +# lambda_n, x_n = var(wm, n, :lambda_n_pipe), var(wm, n, :x_n_pipe) +# +# # Get the corresponding negative flow partitioning (negated). +# partition_n = sort(-get_pipe_flow_partition_negative(ref(wm, n, :pipe, a))) +# bn_range, bn_range_m1 = 1:length(partition_n), 1:length(partition_n)-1 +# +# # Add constraints for the negative flow piecewise approximation. +# c_9 = JuMP.@constraint(wm.model, sum(lambda_n[a, k] for k in bn_range) == 1.0 - y) +# qn_sum = sum(partition_n[k] * lambda_n[a, k] for k in bn_range) +# scalar = _get_scaling_factor(vcat(qn_sum.terms.vals, [1.0])) +# c_10 = JuMP.@constraint(wm.model, scalar * qn_sum == scalar * qn) +# append!(con(wm, n, :pipe_head_loss, a), [c_9, c_10]) +# +# if length(partition_n) > 1 +# # If there are multiple points, constrain the convex combination. +# c_11 = JuMP.@constraint(wm.model, sum(x_n[a, k] for k in bn_range_m1) == 1.0 - y) +# c_12 = JuMP.@constraint(wm.model, lambda_n[a, 1] <= x_n[a, 1]) +# c_13 = JuMP.@constraint(wm.model, lambda_n[a, bn_range[end]] <= x_n[a, bn_range_m1[end]]) +# append!(con(wm, n, :pipe_head_loss, a), [c_11, c_12, c_13]) +# end +# +# # Add a constraint that upper-bounds the head loss variable. +# if maximum(partition_n) != 0.0 +# f_n = r .* partition_n.^exponent +# f_n_ub_expr = sum(f_n[k] * lambda_n[a, k] for k in bn_range) +# scalar = _get_scaling_factor(vcat(f_n_ub_expr.terms.vals, [1.0 / L])) +# c_14 = JuMP.@constraint(wm.model, scalar * dhn / L <= scalar * f_n_ub_expr) +# append!(con(wm, n, :pipe_head_loss, a), [c_14]) +# else +# c_14 = JuMP.@constraint(wm.model, dhn == 0.0) +# append!(con(wm, n, :pipe_head_loss, a), [c_14]) +# end +# +# # Loop over consequential points (i.e., those that have nonzero head loss). +# for flow_value in filter(x -> x > 0.0, partition_n) +# # Add head loss outer (i.e., lower) approximations. +# lhs_n = r * _calc_head_loss_oa(qn, 1.0 - y, flow_value, exponent) +# +# if minimum(abs.(lhs_n.terms.vals)) >= 1.0e-4 +# scalar = _get_scaling_factor(vcat(lhs_n.terms.vals, [1.0 / L])) +# c_15 = JuMP.@constraint(wm.model, scalar * lhs_n <= scalar * dhn / L) +# append!(con(wm, n, :pipe_head_loss, a), [c_15]) +# end +# end +# +# for k in 2:length(partition_n)-1 +# # Add the adjacency constraints for piecewise variables. +# c_16 = JuMP.@constraint(wm.model, lambda_n[a, k] <= x_n[a, k-1] + x_n[a, k]) +# append!(con(wm, n, :pipe_head_loss, a), [c_16]) +# end +# end diff --git a/src/io/common.jl b/src/io/common.jl index ad89bf43..d7bf3e90 100644 --- a/src/io/common.jl +++ b/src/io/common.jl @@ -12,9 +12,9 @@ dictionary of data). Here, `skip_correct` will skip data correction routines (e.g., component status propagation) if set to `true`, and `per_unit` will translate the data model to a per-unit measurement system if set to `true`. """ -function parse_file(path::String; skip_correct::Bool = false, per_unit::Bool = true) +function parse_file(path::String; ne_path::String = "", skip_correct::Bool = false, per_unit::Bool = true) if endswith(path, ".inp") - network_data = WaterModels.parse_epanet(path) + network_data = WaterModels.parse_epanet(path,ne_filename = ne_path) elseif endswith(path, ".json") network_data = WaterModels.parse_json(path) correct_enums!(network_data) @@ -70,7 +70,9 @@ function correct_network_data!(data::Dict{String, <:Any}; per_unit::Bool = true) correct_pipes!(data) correct_des_pipes!(data) correct_pumps!(data) - correct_ne_pumps!(data) + if(haskey(data,"ne_pump")) + correct_ne_pumps!(data) + end correct_regulators!(data) correct_short_pipes!(data) correct_ne_short_pipes!(data) diff --git a/src/io/epanet.jl b/src/io/epanet.jl index ed7733bb..3adeb7a9 100644 --- a/src/io/epanet.jl +++ b/src/io/epanet.jl @@ -4,6 +4,13 @@ function _initialize_epanet_dictionary() return data # Return an empty dictionary of the sectional data. end +function _initialize_epanet_ne_dictionary(data::Dict{String,<:Any}) + ne_data = Dict{String,Any}("time_series" => Dict{String,Any}(), "top_comments" => []) + ne_data["section"] = Dict{String,Array}(x => [] for x in _NE_INP_SECTIONS) + merge!(data["section"],ne_data["section"]) + return data # Return an empty dictionary of the sectional data. +end + function _read_epanet_sections(file_path::String) # Initialize the dictionary that stores raw EPANET data. @@ -54,6 +61,55 @@ function _read_epanet_sections(file_path::String) return data # Return the "raw" EPANET data dictionary. end +function _read_epanet_ne_sections(file_path::String, data::Dict{String,<:Any}) + # Initialize the dictionary that stores raw EPANET NE data. + data = _initialize_epanet_ne_dictionary(data) + + # Instantiate metadata used in parsing EPANET NE sections. + section, line_number = nothing, 0 + + # Populate sections in the NE EPANET dictionary. + for line in eachline(file_path) + # Update important data associated with the line. + line, line_number = strip(line), line_number + 1 + + if length(line) == 0 + # If the line is empty, continue to the next line. + continue + elseif startswith(line, "[") + # If the line starts with "[", it must be a section. + values = split(line; limit = 1) + section_tmp = uppercase(values[1]) + + if section_tmp in _NE_INP_SECTIONS + # If the section title is valid, store it as the current. + section = section_tmp + continue + elseif section_tmp == "[END]" + # If the section title is "[END]", parsing is complete. + break + else + # If the section title is not valid, throw an error. + msg = "File \"$(file_path)\", line $(line_number) has invalid section \"$(section_tmp)\"." + Memento.error(_LOGGER, msg) + end + elseif section === nothing && startswith(line, ";") + # If there is no current section and the section is a comment, store it. + data["top_comments"] = vcat(data["top_comments"], line[2:end]) + continue + elseif section === nothing && !startswith(line, ";") + # If there is no current section and the line is not a comment, throw an error. + msg = "File \"$(file_path)\", line $(line_number) has invalid syntax." + Memento.error(_LOGGER, msg) + end + + # Append a tuple of the line number and the line string to the dictionary. + push!(data["section"][section], (line_number, line)) + end + + return data # Return the "raw" EPANET NE data dictionary. +end + """ parse_epanet(path::String) @@ -68,10 +124,12 @@ original EPANET model, e.g., each pipe with a valve (check or shutoff) is transformed into a WaterModels pipe component _and_ a valve component. Does not perform data correction nor per-unit translations of the data model. """ -function parse_epanet(filename::String) +function parse_epanet(filename::String; ne_filename::String = "") # Read in raw sectional EPANET data. epanet_data = _read_epanet_sections(filename) - + if(!isempty(ne_filename)) + epanet_data = _read_epanet_ne_sections(ne_filename,epanet_data) + end # Parse [OPTIONS] section. _parse_epanet_options(epanet_data) @@ -93,6 +151,11 @@ function parse_epanet(filename::String) # Parse [TANKS] section. _read_tank!(epanet_data) + # Parse [NE_TANKS] section. + if(!isempty(ne_filename)) + _read_ne_tank!(epanet_data) + end + # Add nodes for components with nodal properties. _add_nodes!(epanet_data) @@ -120,9 +183,18 @@ function parse_epanet(filename::String) # Set up a dictionary containing design pipe objects. epanet_data["des_pipe"] = Dict{String,Any}() + # Parse [NE_SHORT_PIPES] section. + if(!isempty(ne_filename)) + _read_ne_short_pipe!(epanet_data) + end + # Parse [PUMPS] section. _read_pump!(epanet_data) + # Parse [NE_PUMPS] section. + if(!isempty(ne_filename)) + _read_ne_pump!(epanet_data) + end # Parse [VALVES] section. _read_regulator!(epanet_data) @@ -138,6 +210,11 @@ function parse_epanet(filename::String) # Update pump data using data from the [ENERGY] section. _update_pump_energy!(epanet_data) + # Update expansion pump data using data from the [ENERGY] section. + if(!isempty(ne_filename)) + _update_ne_pump_energy!(epanet_data) + end + # Parse [COORDINATES] section. _read_coordinate!(epanet_data) @@ -166,6 +243,7 @@ function parse_epanet(filename::String) "duration", "start_time", "demand_charge", + "ne_demand_charge" ] key in keys(epanet_data) && delete!(epanet_data, key) end @@ -494,6 +572,7 @@ function _read_curve!(data::Dict{String,<:Any}) end end + function _read_energy!(data::Dict{String,<:Any}) # Loop over all lines in the [ENERGY] section and parse each. for (line_number, line) in data["section"]["[ENERGY]"] @@ -553,12 +632,54 @@ function _read_energy!(data::Dict{String,<:Any}) else Memento.warn(_LOGGER, "Unknown entry in ENERGY section: $(line)") end + + elseif uppercase(current[1]) == "NE_PUMP" + # Get the corresponding pump data object. + ne_pump_key = findfirst(x -> current[2] == x["source_id"][2], data["ne_pump"]) + ne_pump = data["ne_pump"][ne_pump_key] + + if uppercase(current[3]) == "PRICE" + # Parse the cost of pump operation. + price = parse(Float64, current[4]) # Price per kilowatt hour. + ne_pump["energy_price"] = price * inv(3.6e6) # Price per Joule. + elseif uppercase(current[3]) == "PATTERN" + # Read in the pattern for scaling the pump's energy price. + ne_pump["energy_pattern"] = current[4] + elseif uppercase(current[3]) in ["EFFIC", "EFFICIENCY"] + # Obtain and scale head-versus-flow pump curve. + x = first.(data["curve"][current[4]]) # Flow rate. + y = 0.01 .* last.(data["curve"][current[4]]) # Efficiency. + + if data["flow_units"] == "LPS" # If liters per second... + # Convert from liters per second to cubic meters per second. + x *= 1.0e-3 + elseif data["flow_units"] == "CMH" # If cubic meters per hour... + # Convert from cubic meters per hour to cubic meters per second. + x *= inv(3600.0) + elseif data["flow_units"] == "GPM" # If gallons per minute... + # Convert from gallons per minute to cubic meters per second. + x *= 6.30902e-5 + else + Memento.error(_LOGGER, "Could not find a valid \"units\" option type.") + end + + for k = 1:length(y) + y[k] = x[k] > 0.0 && y[k] > 1.0e-2 ? y[k] : 1.0e-2 + end + + # Curve of efficiency (unitless) versus the flow rate through a pump. + ne_pump["efficiency_curve"] = Array([(x[j], y[j]) for j = 1:length(x)]) + else + Memento.warn(_LOGGER, "Unknown entry in ENERGY section: $(line)") + end else Memento.warn(_LOGGER, "Unknown entry in ENERGY section: $(line)") end end end + + function _read_status!(data::Dict{String,<:Any}) # Loop over all lines in the [STATUS] section and parse each. for (line_number, line) in data["section"]["[STATUS]"] @@ -634,6 +755,54 @@ function _update_pump_energy!(data::Dict{String,<:Any}) end end +function _update_ne_pump_energy!(data::Dict{String,<:Any}) + # Use data previously parsed to infer missing pump data. + for (ne_pump_id, ne_pump) in data["ne_pump"] + # If the pump does not have an efficiency curve, assume a default value. + if !("efficiency_curve" in keys(ne_pump)) && "energy_efficiency" in keys(data) + # If the pump does not have an efficiency, assume the global option. + ne_pump["efficiency"] = data["energy_efficiency"] + elseif !("efficiency_curve" in keys(ne_pump)) + # Otherwise, assume the pump is perfectly efficient. + Memento.warn( + _LOGGER, + "Efficiency for ne pump \"$(ne_pump["name"])\" could not be found.", + ) + ne_pump["efficiency"] = 1.0 + end + + # If the pump does not have an energy pattern, assume the global option. + if !("energy_pattern" in keys(ne_pump)) && "energy_pattern" in keys(data) + ne_pump["energy_pattern"] = data["energy_pattern"] + end + + # If the energy price is not specified for the pump, set a reasonable value. + if !("energy_price" in keys(ne_pump)) && "energy_price" in keys(data) + # If the pump does not have an energy price, assume the global option. + ne_pump["energy_price"] = data["energy_price"] + elseif !("energy_price" in keys(ne_pump)) # If there is no energy price... + # Otherwise, assume the pump does not have an energy price. + Memento.warn(_LOGGER, "Price for ne pump \"$(ne_pump["name"])\" could not be found.") + ne_pump["energy_price"] = 0.0 + end + + # If an energy pattern is specified, store the time series of prices. + if "energy_pattern" in keys(ne_pump) + # Build the pattern of prices using existing data. + pattern_name = ne_pump["energy_pattern"] + price_pattern = ne_pump["energy_price"] .* data["pattern"][pattern_name] + push!(price_pattern, minimum(price_pattern)) + + # Add the varying price data to the pump time series entry. + entry = Dict{String,Array{Float64}}("energy_price" => price_pattern) + data["time_series"]["ne_pump"][ne_pump_id] = entry + + # Delete the "energy_pattern" data from the pump. + delete!(ne_pump, "energy_pattern") + end + end +end + "Parse and store demand data from an EPANET-formatted data dictionary." function _read_demand!(data::Dict{String,<:Any}) # Initialize dictionaries associated with demands. @@ -731,7 +900,6 @@ function _add_node_map!(data::Dict{String,<:Any}) map = Dict{String,Int}(x["name"] => x["node"] for (i, x) in data[type]) data["node_map"] = merge(data["node_map"], map) end - # Ensure the number of node-type components is equal to the number of nodes. num_components = sum(length(data[type]) for type in ["demand", "reservoir", "tank"]) @@ -906,6 +1074,36 @@ function _read_pipe!(data::Dict{String,<:Any}) data["pipe"] = _transform_component_indices(data["pipe"]) end +function _read_ne_short_pipe!(data::Dict{String,<:Any}) + # Initialize dictionary associated with expansion short pipes. + data["ne_short_pipe"] = Dict{String,Dict{String,Any}}() + + # Initialize a temporary index to be updated while parsing. + index::Int = 0 + + # Loop over all lines in the [PIPES] section and parse each. + for (line_number, line) in data["section"]["[NE_SHORT_PIPES]"] + current = split(split(line, ";")[1]) + length(current) == 0 && continue # Skip. + + # Initialize the pipe entry to be added. + ne_short_pipe = Dict{String,Any}("name" => current[1], "status" => STATUS_UNKNOWN) + ne_short_pipe["source_id"] = ["ne_short_pipe", current[1]] + ne_short_pipe["node_fr"] = data["node_map"][current[2]] + ne_short_pipe["node_to"] = data["node_map"][current[3]] + ne_short_pipe["construction_cost"] = parse(Float64, current[5]) + + # Add a temporary index to be used in the data dictionary. + ne_short_pipe["index"] = string(index += 1) + + # Append the pipe entry to the data dictionary. + data["ne_short_pipe"][current[1]] = ne_short_pipe + end + + # Replace with new dictionaries that use integer component indices. + data["ne_short_pipe"] = _transform_component_indices(data["ne_short_pipe"]) +end + function _read_pump!(data::Dict{String,<:Any}) # Initialize dictionaries associated with pumps. data["pump"] = Dict{String,Dict{String,Any}}() @@ -963,7 +1161,7 @@ function _read_pump!(data::Dict{String,<:Any}) end # Curve of head (meters) versus flow (cubic meters per second). - pump["head_curve"] = Array([(flow[j], head[j]) for j = 1:length(flow)]) + 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 @@ -985,6 +1183,87 @@ function _read_pump!(data::Dict{String,<:Any}) data["pump"], data["time_series"]["pump"]) end +function _read_ne_pump!(data::Dict{String,<:Any}) + # Initialize dictionaries associated with expansion pumps. + data["ne_pump"] = Dict{String,Dict{String,Any}}() + data["time_series"]["ne_pump"] = Dict{String,Any}() + + # Initialize a temporary index to be updated while parsing. + index::Int = 0 + + # Loop over all lines in the [NE_PUMPS] section and parse each. + for (line_number, line) in data["section"]["[NE_PUMPS]"] + current = split(split(line, ";")[1]) + length(current) == 0 && continue # Skip. + + # Initialize the pipe entry to be added. + ne_pump = Dict{String,Any}("name" => current[1], "status" => STATUS_UNKNOWN) + ne_pump["source_id"] = ["ne_pump", current[1]] + ne_pump["node_fr"] = data["node_map"][current[2]] + ne_pump["node_to"] = data["node_map"][current[3]] + ne_pump["construction_cost"] = parse(Float64,current[length(current)]) + + # Loop over remaining entries and store remaining properties. + for i in range(4; stop = length(current)-1, step = 2) + if uppercase(current[i]) == "POWER" + # Memento.error(_LOGGER, "Constant power pumps are not supported.") + ne_pump["pump_type"] = PUMP_CONSTANT_POWER + + # Parse the head of the reservoir (in meters). + if data["flow_units"] == "LPS" || data["flow_units"] == "CMH" + # Conver power from kilowatts to watts. + ne_pump["power"] = 1.0e3 * parse(Float64, current[i+1]) + elseif data["flow_units"] == "GPM" # If gallons per minute... + # Convert power from horsepower to watts. + ne_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 + + # 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). + ne_pump["head_curve"] = Array([(flow[j], head[j]) for j = 1:length(flow)]) + elseif uppercase(current[i]) != "HEAD" + Memento.error(_LOGGER, "Ne Pump keyword in INP file not recognized.") + end + end + + # Flow is always in the positive direction for pumps. + ne_pump["flow_direction"] = FLOW_DIRECTION_POSITIVE + + # Add a temporary index to be used in the data dictionary. + ne_pump["index"] = string(index += 1) + + # Append the ne_pump entry to the data dictionary. + data["ne_pump"][current[1]] = ne_pump + end + + # Replace with new dictionaries that use integer component indices. + data["ne_pump"] = _transform_component_indices(data["ne_pump"]) + data["time_series"]["ne_pump"] = _transform_time_series_indices( + data["ne_pump"], data["time_series"]["ne_pump"]) +end + + function _read_reservoir!(data::Dict{String,<:Any}) # Initialize dictionaries associated with reservoirs. data["reservoir"] = Dict{String,Dict{String,Any}}() @@ -1112,6 +1391,62 @@ function _read_tank!(data::Dict{String,<:Any}) data["tank"] = _transform_component_indices(data["tank"]) end +function _read_ne_tank!(data::Dict{String,<:Any}) + # Initialize dictionaries associated with tanks. + data["ne_tank"] = Dict{String,Dict{String,Any}}() + + # Initialize a temporary index to be updated while parsing. + index::Int = length(data["tank"]) + + # Loop over all lines in the [TANKS] section and parse each. + for (line_number, line) in data["section"]["[NE_TANKS]"] + current = split(split(line, ";")[1]) + length(current) == 0 && continue # Skip. + + # Initialize the tank entry to be added. + ne_tank = Dict{String,Any}("name" => current[1], "status" => STATUS_ACTIVE) + ne_tank["source_id"] = ["tank", current[1]] + ne_tank["dispatchable"] = false + + # Store all measurements associated with expansion tanks in metric units. + if data["flow_units"] == "LPS" || data["flow_units"] == "CMH" + # Retain the original length values (in meters). + ne_tank["elevation"] = parse(Float64, current[2]) + ne_tank["init_level"] = parse(Float64, current[3]) + ne_tank["min_level"] = parse(Float64, current[4]) + ne_tank["max_level"] = parse(Float64, current[5]) + ne_tank["diameter"] = parse(Float64, current[6]) + + # Retain the original minimum volume (in cubic meters). + ne_tank["min_vol"] = parse(Float64, current[7]) + elseif data["flow_units"] == "GPM" # If gallons per minute... + # Convert length values from feet to meters. + ne_tank["elevation"] = 0.3048 * parse(Float64, current[2]) + ne_tank["init_level"] = 0.3048 * parse(Float64, current[3]) + ne_tank["min_level"] = 0.3048 * parse(Float64, current[4]) + ne_tank["max_level"] = 0.3048 * parse(Float64, current[5]) + ne_tank["diameter"] = 0.3048 * parse(Float64, current[6]) + + # Convert minimum volume from cubic feet to cubic meters. + ne_tank["min_vol"] = 0.3048^3 * parse(Float64, current[7]) + else + Memento.error(_LOGGER, "Could not find a valid \"Units\" option type.") + end + + # Add a temporary index to be used in the data dictionary. + ne_tank["index"] = string(index += 1) + + # Append the tank entry to the data dictionary. + data["ne_tank"][current[1]] = ne_tank + end + + # Replace with new dictionaries that use integer component indices. + data["ne_tank"] = _transform_component_indices(data["ne_tank"]) + println("Warning: Temporarily labeled the source_id of ne_tanks as tanks and merged them into regular tanks") + merge!(data["tank"],data["ne_tank"]) + delete!(data, "ne_tank") +end + function _read_time!(data::Dict{String,<:Any}) # Loop over all lines in the [TIMES] section and parse each. for (line_number, line) in data["section"]["[TIMES]"] @@ -1402,3 +1737,9 @@ _INP_SECTIONS = [ "[BACKDROP]", "[TAGS]", ] + +_NE_INP_SECTIONS = [ + "[NE_PUMPS]", + "[NE_TANKS]", + "[NE_SHORT_PIPES]", +] diff --git a/src/prob/ne.jl b/src/prob/ne.jl index a57097ac..ad4c9bd5 100644 --- a/src/prob/ne.jl +++ b/src/prob/ne.jl @@ -27,10 +27,18 @@ function build_mn_ne(wm::AbstractWaterModel) # Ensure tanks recover their initial volume. n_1, n_f = network_ids[1], network_ids[end] - for i in ids(wm, n_f, :tank) - constraint_tank_volume_recovery(wm, i, n_1, n_f) + if(haskey(wm.ref[:it][wm_it_sym],:tank_volume_recovery_time_points)) + tank_volume_recovery_time_points = wm.ref[:it][wm_it_sym][:tank_volume_recovery_time_points] + else + tank_volume_recovery_time_points = Set([]) + end + union(tank_volume_recovery_time_points ,n_f) + for n_tank in tank_volume_recovery_time_points + for i in ids(wm, n_tank, :tank) + constraint_tank_volume_recovery(wm, i, n_1, n_tank) + end end # Add the optimal network expansion objective. objective_ne(wm) -end \ No newline at end of file +end diff --git a/src/prob/owf.jl b/src/prob/owf.jl index ab7b5052..eb59ec22 100644 --- a/src/prob/owf.jl +++ b/src/prob/owf.jl @@ -32,8 +32,16 @@ function build_mn_owf(wm::AbstractWaterModel) # Ensure tanks recover their initial volume. n_1, n_f = network_ids[1], network_ids[end] - for i in ids(wm, n_f, :tank) - constraint_tank_volume_recovery(wm, i, n_1, n_f) + if(haskey(wm.ref[:it][wm_it_sym],:tank_volume_recovery_time_points)) + tank_volume_recovery_time_points = wm.ref[:it][wm_it_sym][:tank_volume_recovery_time_points] + else + tank_volume_recovery_time_points = Set([]) + end + union(tank_volume_recovery_time_points ,n_f) + for n_tank in tank_volume_recovery_time_points + for i in ids(wm, n_tank, :tank) + constraint_tank_volume_recovery(wm, i, n_1, n_tank) + end end # Add the optimal water flow objective. diff --git a/src/prob/wf.jl b/src/prob/wf.jl index 5d89526a..b2efbe63 100644 --- a/src/prob/wf.jl +++ b/src/prob/wf.jl @@ -15,20 +15,24 @@ end function build_wf(wm::AbstractWaterModel) # Create head loss functions, if necessary. - _function_head_loss(wm) + # _function_head_loss(wm) # Physical variables. variable_head(wm) variable_flow(wm) variable_pump_head_gain(wm) + variable_ne_pump_head_gain(wm) variable_pump_power(wm) + variable_ne_pump_power(wm) # Indicator (status) variables. variable_des_pipe_indicator(wm) variable_pump_indicator(wm) + variable_ne_pump_indicator(wm) variable_regulator_indicator(wm) variable_valve_indicator(wm) variable_ne_short_pipe_indicator(wm) + variable_ne_pump_build(wm) # Create flow-related variables for node attachments. variable_demand_flow(wm) @@ -67,6 +71,7 @@ function build_wf(wm::AbstractWaterModel) constraint_on_off_pump_head(wm, a) constraint_on_off_pump_head_gain(wm, a) constraint_on_off_pump_flow(wm, a) + # println("*******Disabling pump power constraint*****************") constraint_on_off_pump_power(wm, a) end @@ -76,6 +81,7 @@ function build_wf(wm::AbstractWaterModel) constraint_on_off_pump_head_gain_ne(wm, a) constraint_on_off_pump_flow_ne(wm, a) constraint_on_off_pump_power_ne(wm, a) + constraint_on_off_pump_build_ne(wm, a) end for (k, pump_group) in ref(wm, :pump_group) @@ -119,7 +125,7 @@ end function build_mn_wf(wm::AbstractWaterModel) # Create head loss functions, if necessary. - _function_head_loss(wm) + # _function_head_loss(wm) # Get all network IDs in the multinetwork. network_ids = sort(collect(nw_ids(wm))) @@ -129,7 +135,7 @@ function build_mn_wf(wm::AbstractWaterModel) else network_ids_inner = network_ids end - + # println("Running SCIP friendly (LP) version ******") for n in network_ids_inner # Physical variables. variable_head(wm; nw=n) @@ -140,7 +146,6 @@ function build_mn_wf(wm::AbstractWaterModel) variable_ne_pump_power(wm;nw=n) - # Indicator (status) variables. variable_des_pipe_indicator(wm; nw=n) variable_pump_indicator(wm; nw=n) @@ -155,6 +160,15 @@ function build_mn_wf(wm::AbstractWaterModel) variable_reservoir_flow(wm; nw=n) variable_tank_flow(wm; nw=n) + # + if(haskey(wm.ref[:it][wm_it_sym],:sol_from_relaxation)) + sol_from_relaxation = wm.ref[:it][wm_it_sym][:sol_from_relaxation] + if(n == 1) + println("Fixing Binary Values from relaxed solution ----------------------------") + end + fix_variables_to_relaxed_solutions(wm, sol_from_relaxation; nw=n) + end + # Flow conservation at all nodes. for i in ids(wm, :node; nw=n) constraint_flow_conservation(wm, i; nw=n) @@ -166,8 +180,8 @@ function build_mn_wf(wm::AbstractWaterModel) constraint_pipe_flow(wm, a; nw=n) constraint_pipe_head(wm, a; nw=n) constraint_pipe_head_loss(wm, a; nw=n) + # constraint_pipe_head_loss_scip_lp(wm, a; nw=n) end - # Ensure design pipes are not present in the problem. @assert length(ids(wm, :des_pipe_arc; nw=n)) == 0 @assert length(ids(wm, :des_pipe; nw=n)) == 0 @@ -222,12 +236,12 @@ function build_mn_wf(wm::AbstractWaterModel) constraint_on_off_valve_head(wm, a; nw=n) constraint_on_off_valve_flow(wm, a; nw=n) end + end # Start with the first network, representing the initial time step. n_1 = network_ids[1] - # Constraints on tank volumes. for i in ids(wm, :tank; nw = n_1) # Set initial conditions of tanks. constraint_tank_volume(wm, i; nw = n_1) @@ -264,6 +278,7 @@ function build_mn_wf(wm::AbstractWaterModel) end end + # println(wm.model) # Add the objective. objective_wf(wm) end diff --git a/src/util/fixing_variables.jl b/src/util/fixing_variables.jl new file mode 100644 index 00000000..e596b37d --- /dev/null +++ b/src/util/fixing_variables.jl @@ -0,0 +1,168 @@ +function fix_variables_to_relaxed_solutions(wm::AbstractWaterModel, sol_val_dict::Dict{String, Any}; nw::Int=nw_id_default) + + + + + fixing_vector = sol_val_dict["fixing_vector"] + fixing_all_pipes_direction = fixing_vector[1] + fixing_pumps_direction = fixing_vector[2] + fixing_build = fixing_vector[3] + fixing_activation = fixing_vector[4] + + understanding_partitions = true + + if(fixing_all_pipes_direction == true) + if(haskey(ref(wm,nw),:pipe)) + if(nw == 1) + println("fixing pipe flow direction") + end + for (a, pipe) in ref(wm, nw, :pipe) + val = sol_val_dict["nw"][string(nw)]["pipe"][string(a)]["y"] + # println("fixing pipe flow direction nw $nw pipe $a to $val") + y_pipe = var(wm, nw, :y_pipe, a) + JuMP.unset_binary(y_pipe) + JuMP.fix(y_pipe, val; force = true) + + + # if(understanding_partitions == true) + # lambda_pipe = var(wm, nw, :lambda_p_pipe) + # x_pipe = var(wm, nw, :x_p_pipe) + # end + end + if(haskey(ref(wm,nw),:short_pipe)) + if(nw == 1) + println("fixing short pipe flow direction") + end + for (a, short_pipe) in ref(wm, nw, :short_pipe) + val = sol_val_dict["nw"][string(nw)]["short_pipe"][string(a)]["y"] + val = Int(round(val)) + y_short_pipe = var(wm, nw, :y_short_pipe, a) + JuMP.unset_binary(y_short_pipe) + JuMP.fix(y_short_pipe, val; force = true) + end + end + + if(haskey(ref(wm,nw),:ne_short_pipe)) + if(nw == 1) + println("fixing ne short pipe direction") + end + for (a, ne_short_pipe) in ref(wm, nw, :ne_short_pipe) + val = sol_val_dict["nw"][string(nw)]["ne_short_pipe"][string(a)]["y"] + val = Int(round(val)) + y_ne_short_pipe = var(wm, nw, :y_ne_short_pipe, a) + JuMP.unset_binary(y_ne_short_pipe) + JuMP.fix(y_ne_short_pipe, val; force = true) + end + end + + end + end + + + if(fixing_pumps_direction == true) + + if(haskey(ref(wm,nw),:pump)) + if(nw == 1) + println("fixing pump flow direction") + end + for (a, pump) in ref(wm, nw, :pump) + val = sol_val_dict["nw"][string(nw)]["pump"][string(a)]["y"] + val = Int(round(val)) + y_pump = var(wm, nw, :y_pump, a) + JuMP.unset_binary(y_pump) + JuMP.fix(y_pump, val; force = true) + end + end + if(haskey(ref(wm,nw),:ne_pump)) + if(nw == 1) + println("fixing ne pump flow direction") + end + for (a, ne_pump) in ref(wm, nw, :ne_pump) + val = sol_val_dict["nw"][string(nw)]["ne_pump"][string(a)]["y"] + val = Int(round(val)) + y_ne_pump = var(wm, nw, :y_ne_pump, a) + JuMP.unset_binary(y_ne_pump) + JuMP.fix(y_ne_pump, val; force = true) + end + end + + # if(haskey(ref(wm,nw),:valve) && haskey(sol_val_dict["nw"][string(nw)], "valve")) + # for (a, valve) in ref(wm, nw, :valve) + # yval = sol_val_dict["nw"][string(nw)]["valve"][string(a)]["y"] + # yval = Int(round(yval)) + # y_valve = var(wm, nw, :y_valve, a) + # JuMP.unset_binary(y_valve) + # JuMP.fix(y_valve, yval; force = true) + # end + # end + end + + if(fixing_build == true) + if(haskey(ref(wm,nw),:ne_pump)) + if(nw == 1) + println("fixing ne pump build status") + end + for (a, ne_pump) in ref(wm, nw, :ne_pump) + val = sol_val_dict["nw"][string(nw)]["ne_pump"][string(a)]["build_status"] + val = Int(round(val)) + x_ne_pump = var(wm, nw, :x_ne_pump, a) + JuMP.unset_binary(x_ne_pump) + JuMP.fix(x_ne_pump, val; force = true) + end + end + if(haskey(ref(wm,nw),:ne_short_pipe)) + if(nw == 1) + println("fixing ne short pipe build status") + end + for (a, ne_short_pipe) in ref(wm, nw, :ne_short_pipe) + val = sol_val_dict["nw"][string(nw)]["ne_short_pipe"][string(a)]["status"] + val = Int(round(val)) + z_ne_short_pipe = var(wm, nw, :z_ne_short_pipe, a) + JuMP.unset_binary(z_ne_short_pipe) + JuMP.fix(z_ne_short_pipe, val; force = true) + end + end + end + + if(fixing_activation == true) + if(haskey(ref(wm,nw),:pump)) + if(nw == 1) + println("fixing pump activation status") + end + for (a, pump) in ref(wm, nw, :pump) + val = sol_val_dict["nw"][string(nw)]["pump"][string(a)]["status"] + val = Int(round(val)) + z_pump = var(wm, nw, :z_pump, a) + JuMP.unset_binary(z_pump) + JuMP.fix(z_pump, val; force = true) + end + end + + if(haskey(ref(wm,nw),:ne_pump)) + if(nw == 1) + println("fixing ne pump activation status") + end + for (a, ne_pump) in ref(wm, nw, :ne_pump) + val = sol_val_dict["nw"][string(nw)]["ne_pump"][string(a)]["status"] + val = Int(round(val)) + z_ne_pump = var(wm, nw, :z_ne_pump, a) + JuMP.unset_binary(z_ne_pump) + JuMP.fix(z_ne_pump, val; force = true) + end + end + end + + + # if(haskey(ref(wm,nw),:valve) && haskey(sol_val_dict["nw"][string(nw)], "valve")) + # for (a, valve) in ref(wm, nw, :valve) + # val = sol_val_dict["nw"][string(nw)]["valve"][string(a)]["status"] + # val = Int(round(val)) + # z_valve = var(wm, nw, :z_valve, a) + # JuMP.unset_binary(z_valve) + # JuMP.fix(z_valve, val; force = true) + # end + + # end + + +end diff --git a/test/data/epanet/snapshot/balerma.inp b/test/data/epanet/snapshot/balerma.inp index 18b32df9..c4a5e4ac 100644 --- a/test/data/epanet/snapshot/balerma.inp +++ b/test/data/epanet/snapshot/balerma.inp @@ -3,9 +3,9 @@ Balerma Network [JUNCTIONS] ;ID Elev Demand Pattern - 179001 60.0000 5.550000 ; - 179 60.0000 5.550000 ; - 177 58.9000 5.550000 ; + 179001 60.0000 5.550000 ; + 179 60.0000 5.550000 sk_dem ; + 177 58.9000 5.550000 sk_dem ; 174 60.0000 5.550000 ; 173 59.3000 5.550000 ; 171001 60.0000 5.550000 ; @@ -19,7 +19,7 @@ Balerma Network 163 55.2000 5.550000 ; 164 57.8000 5.550000 ; 162 55.2000 5.550000 ; - 161 55.0000 5.550000 ; + 161 55.0000 5.550000 sk_dem ; 106 54.0000 5.550000 ; 124 52.2000 5.550000 ; 125 51.1000 5.550000 ; @@ -930,6 +930,7 @@ Balerma Network [PATTERNS] ;ID Multipliers +sk_dem 1 1 1 1 1 1 [CURVES] ;ID X-Value Y-Value @@ -973,7 +974,7 @@ Balerma Network ;Tank Model [TIMES] - Duration 0:00 + Duration 6:00 Hydraulic Timestep 1:00 Quality Timestep 0:06 Pattern Timestep 1:00