diff --git a/CHANGELOG.md b/CHANGELOG.md index 798a818f..cae881c0 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,6 +1,14 @@ WaterModels.jl Change Log ========================= +### v0.6.0 +- Decoupled `check_valve` and `shutoff_valve` objects from `pipe` objects. +- Introduced a new `valve` component. +- Introduced a new `short_pipe` component. +- Renamed `pressure_reducing_valve` to `regulator`. +- Decomposed component constraints a bit more. +- Implemented `QRD` and `CQRD` formulations. + ### v0.5.0 - Rename formulation types. diff --git a/README.md b/README.md index 53948ee2..1f2d88a8 100644 --- a/README.md +++ b/README.md @@ -18,6 +18,8 @@ This decoupling enables the definition of a wide variety of optimization formula * CRD - continuous (convex) relaxation- and direction-based formulation * LA - linear approximation-based formulation * LRD - linear relaxation- and direction-based formulation +* QRD - nonconvex quadratic relaxation- and direction-based formulation +* CQRD - convex quadratic relaxation- and direction-based formulation ## Documentation The package [documentation](https://lanl-ansi.github.io/WaterModels.jl/latest) includes a [quick start guide](https://lanl-ansi.github.io/WaterModels.jl/latest/quickguide). diff --git a/src/WaterModels.jl b/src/WaterModels.jl index 93072601..b02ce4a2 100644 --- a/src/WaterModels.jl +++ b/src/WaterModels.jl @@ -61,6 +61,8 @@ include("form/crd.jl") include("form/la.jl") include("form/lrd.jl") include("form/nc.jl") +include("form/qrd.jl") +include("form/cqrd.jl") include("prob/wf.jl") include("prob/owf.jl") @@ -68,6 +70,7 @@ include("prob/des.jl") include("util/unbinarize.jl") include("util/obbt.jl") +include("util/compute_cuts.jl") include("core/export.jl") end diff --git a/src/core/base.jl b/src/core/base.jl index f9a73988..7b76af21 100644 --- a/src/core/base.jl +++ b/src/core/base.jl @@ -137,13 +137,13 @@ system-wide values that need to be computed globally. Some of the common keys include: -* `:pipe` -- the set of pipes in the network without valves, -* `:check_valve` -- the set of pipes in the network with check valves, -* `:shutoff_valve` -- the set of pipes in the network wiith shutoff valves, -* `:pressure_reducing_valve` -- the set of pressure reducing valves in the network, +* `:pipe` -- the set of pipes in the network, * `:pump` -- the set of pumps in the network, +* `:regulator` -- the set of regulators in the network, +* `:short_pipe` -- the set of short pipes in the network, +* `:valve` -- the set of valves in the network, * `:node` -- the set of all nodes in the network, -* `:junction` -- the set of junctions in the network, +* `:demand` -- the set of demands in the network, * `:reservoir` -- the set of reservoirs in the network, * `:tank` -- the set of tanks in the network """ @@ -155,32 +155,25 @@ end function _ref_add_core!(nw_refs::Dict{Int,<:Any}) for (nw, ref) in nw_refs # Collect dispatchable and nondispatchable nodal components in the network. - ref[:dispatchable_junction] = filter(x -> x.second["dispatchable"], ref[:junction]) - ref[:nondispatchable_junction] = filter(x -> !x.second["dispatchable"], ref[:junction]) + ref[:dispatchable_demand] = filter(x -> x.second["dispatchable"], ref[:demand]) + ref[:nondispatchable_demand] = filter(x -> !x.second["dispatchable"], ref[:demand]) # Compute resistances for pipe-type components in the network. ref[:resistance] = calc_resistances(ref[:pipe], ref[:viscosity], ref[:head_loss]) ref[:resistance_cost] = calc_resistance_costs(ref[:pipe], ref[:viscosity], ref[:head_loss]) - - # TODO: Check and shutoff valves should not have pipe properties. - ref[:check_valve] = filter(x -> x.second["has_check_valve"], ref[:pipe]) - ref[:shutoff_valve] = filter(x -> x.second["has_shutoff_valve"], ref[:pipe]) - ref[:pipe] = filter(x -> !x.second["has_check_valve"], ref[:pipe]) - ref[:pipe] = filter(x -> !x.second["has_shutoff_valve"], ref[:pipe]) ref[:des_pipe] = filter(is_des_pipe, ref[:pipe]) ref[:pipe] = filter(!is_des_pipe, ref[:pipe]) # Create mappings of "from" and "to" arcs for link- (i.e., edge-) type components. - for name in - ["check_valve", "shutoff_valve", "pipe", "pressure_reducing_valve", "pump"] + for name in ["pipe", "pump", "regulator", "short_pipe", "valve"] fr_sym, to_sym = Symbol(name * "_fr"), Symbol(name * "_to") - ref[fr_sym] = [(i, c["node_fr"], c["node_to"]) for (i, c) in ref[Symbol(name)]] - ref[to_sym] = [(i, c["node_to"], c["node_fr"]) for (i, c) in ref[Symbol(name)]] + ref[fr_sym] = [(a, c["node_fr"], c["node_to"]) for (a, c) in ref[Symbol(name)]] + ref[to_sym] = [(a, c["node_to"], c["node_fr"]) for (a, c) in ref[Symbol(name)]] end # Set up dictionaries mapping node indices to attached component indices. - for name in ["junction", "tank", "reservoir"] + for name in ["demand", "tank", "reservoir"] name_sym = Symbol("node_" * name) ref[name_sym] = Dict{Int,Array{Int,1}}(i => Int[] for (i, node) in ref[:node]) diff --git a/src/core/constraint.jl b/src/core/constraint.jl index a14451e9..b9626de1 100644 --- a/src/core/constraint.jl +++ b/src/core/constraint.jl @@ -3,84 +3,87 @@ ######################################################################### -""" - constraint_reservoir_head(wm, n, i, head) -""" -function constraint_reservoir_head(wm::AbstractWaterModel, n::Int, i::Int, head::Float64) - h = var(wm, n, :h, i) - c = JuMP.@constraint(wm.model, h == head) - con(wm, n, :reservoir_head)[i] = c -end - - """ constraint_flow_conservation( - wm, n, i, check_valve_fr, check_valve_to, pipe_fr, pipe_to, des_pipe_fr, - des_pipe_to, pump_fr, pump_to, pressure_reducing_valve_fr, - pressure_reducing_valve_to, shutoff_valve_fr, shutoff_valve_to, reservoirs, tanks, - dispatachable_junctions, fixed_demand) + wm, n, i, pipe_fr, pipe_to, des_pipe_fr, des_pipe_to, pump_fr, pump_to, + regulator_fr, regulator_to, short_pipe_fr, short_pipe_to, valve_fr, valve_to, + reservoirs, tanks, dispatachable_demands, fixed_demand) + +Adds a constraint that ensures flow conservation at a node in the network. """ function constraint_flow_conservation( - wm::AbstractWaterModel, n::Int, i::Int, check_valve_fr::Array{Int64,1}, - check_valve_to::Array{Int64,1}, pipe_fr::Array{Int64,1}, pipe_to::Array{Int64,1}, - des_pipe_fr::Array{Int64,1}, des_pipe_to::Array{Int64,1}, pump_fr::Array{Int64,1}, - pump_to::Array{Int64,1}, pressure_reducing_valve_fr::Array{Int64,1}, - pressure_reducing_valve_to::Array{Int64,1}, shutoff_valve_fr::Array{Int64,1}, - shutoff_valve_to::Array{Int64,1}, reservoirs::Array{Int64,1}, tanks::Array{Int64,1}, - dispatchable_junctions::Array{Int64,1}, fixed_demand::Float64) + wm::AbstractWaterModel, n::Int, i::Int, pipe_fr::Array{Int64,1}, + pipe_to::Array{Int64,1}, des_pipe_fr::Array{Int64,1}, des_pipe_to::Array{Int64,1}, + pump_fr::Array{Int64,1}, pump_to::Array{Int64,1}, regulator_fr::Array{Int64,1}, + regulator_to::Array{Int64,1}, short_pipe_fr::Array{Int64,1}, + short_pipe_to::Array{Int64,1}, valve_fr::Array{Int64,1}, valve_to::Array{Int64,1}, + reservoirs::Array{Int64,1}, tanks::Array{Int64,1}, dispatchable_demands::Array{Int64,1}, + fixed_demand::Float64) # Collect flow variable references per component. - q_check_valve = var(wm, n, :q_check_valve) - q_pipe = var(wm, n, :q_pipe) - q_des_pipe = var(wm, n, :q_des_pipe_sum) - q_pump = var(wm, n, :q_pump) - q_pressure_reducing_valve = var(wm, n, :q_pressure_reducing_valve) - q_shutoff_valve = var(wm, n, :q_shutoff_valve) + q_pipe, q_des_pipe = var(wm, n, :q_pipe), var(wm, n, :q_des_pipe_sum) + q_pump, q_regulator = var(wm, n, :q_pump), var(wm, n, :q_regulator) + q_short_pipe, q_valve = var(wm, n, :q_short_pipe), var(wm, n, :q_valve) + q_reservoir, q_tank = var(wm, n, :q_reservoir), var(wm, n, :q_tank) + q_demand = var(wm, n, :q_demand) # Add the flow conservation constraint. - c = JuMP.@constraint(wm.model, - - sum(q_check_valve[a] for a in check_valve_fr) + - sum(q_check_valve[a] for a in check_valve_to) - + con(wm, n, :flow_conservation)[i] = JuMP.@constraint(wm.model, - sum(q_pipe[a] for a in pipe_fr) + sum(q_pipe[a] for a in pipe_to) - sum(q_des_pipe[a] for a in des_pipe_fr) + sum(q_des_pipe[a] for a in des_pipe_to) - sum(q_pump[a] for a in pump_fr) + sum(q_pump[a] for a in pump_to) - - sum(q_pressure_reducing_valve[a] for a in pressure_reducing_valve_fr) + - sum(q_pressure_reducing_valve[a] for a in pressure_reducing_valve_to) - - sum(q_shutoff_valve[a] for a in shutoff_valve_fr) + - sum(q_shutoff_valve[a] for a in shutoff_valve_to) == - - sum(var(wm, n, :qr, k) for k in reservoirs) - - sum(var(wm, n, :qt, k) for k in tanks) + - sum(var(wm, n, :demand, k) for k in dispatchable_junctions) + fixed_demand) - - con(wm, n, :flow_conservation)[i] = c + sum(q_regulator[a] for a in regulator_fr) + + sum(q_regulator[a] for a in regulator_to) - + sum(q_short_pipe[a] for a in short_pipe_fr) + + sum(q_short_pipe[a] for a in short_pipe_to) - + sum(q_valve[a] for a in valve_fr) + sum(q_valve[a] for a in valve_to) == - + sum(q_reservoir[k] for k in reservoirs) - sum(q_tank[k] for k in tanks) + + sum(q_demand[k] for k in dispatchable_demands) + fixed_demand) end -function constraint_tank_state_initial(wm::AbstractWaterModel, n::Int, i::Int, V_0::Float64) +""" + constraint_tank_volume_fixed(wm, n, i, V_0) + +Adds a constraint that ensures the volume of a tank at some time step is fixed. Here, `wm` +is the WaterModels object, `n` is the index of a subnetwork within a multinetwork, `i` is +the index of the tank, and `V_0` is the fixed volume of the tank that is desired. +""" +function constraint_tank_volume_fixed(wm::AbstractWaterModel, n::Int, i::Int, V_0::Float64) V = var(wm, n, :V, i) c = JuMP.@constraint(wm.model, V == V_0) - con(wm, n, :tank_state)[i] = c + con(wm, n, :tank_volume)[i] = c end """ - constraint_tank_state(wm, n_1, n_2, i, time_step) + constraint_tank_volume(wm, n_1, n_2, i, time_step) Adds a constraint that integrates the volume of a tank forward in time. Here, `wm` is the WaterModels object, `n_1` is the index of a subnetwork within a multinetwork, `n_2` is the index of another subnetwork forward in time, relative to `n_1`, i is the index of the tank, -and time_step is the time step (in seconds) between networks `n_1` and `n_2`. +and time_step is the time step (in seconds) of the interval from network `n_1` to `n_2`. """ -function constraint_tank_state(wm::AbstractWaterModel, n_1::Int, n_2::Int, i::Int, time_step::Float64) - qt = var(wm, n_1, :qt, i) # Tank outflow. +function constraint_tank_volume(wm::AbstractWaterModel, n_1::Int, n_2::Int, i::Int, time_step::Float64) + qt = var(wm, n_1, :q_tank, i) # Tank outflow. V_1, V_2 = var(wm, n_1, :V, i), var(wm, n_2, :V, i) c = JuMP.@constraint(wm.model, V_1 - time_step * qt == V_2) - con(wm, n_2, :tank_state)[i] = c + con(wm, n_2, :tank_volume)[i] = c end -function constraint_recover_volume(wm::AbstractWaterModel, i::Int, n_1::Int, n_f::Int) - _initialize_con_dict(wm, :recover_volume, nw=n_f) - V_1, V_f = var(wm, n_1, :V, i), var(wm, n_f, :V, i) - c = JuMP.@constraint(wm.model, V_f >= V_1) - con(wm, n_f, :recover_volume)[i] = c +""" + constraint_tank_volume_recovery(wm, i, n_1, n_f) + +Adds a constraint that ensures the volume of a tank at the end of the time horizon is +greater than or equal to the volume of the tank at the beginning of the time horizon. Here, +`wm` is the WaterModels object, `n_1` is the index of the first subnetwork within a +multinetwork, `n_f` is the index of the final subnetwork, and i is the index of the tank. +""" +function constraint_tank_volume_recovery(wm::AbstractWaterModel, i::Int, n_1::Int, n_f::Int) + if !ref(wm, n_f, :tank, i)["dispatchable"] + _initialize_con_dict(wm, :tank_volume_recovery, nw=n_f) + V_1, V_f = var(wm, n_1, :V, i), var(wm, n_f, :V, i) + c = JuMP.@constraint(wm.model, V_f >= V_1) + con(wm, n_f, :tank_volume_recovery)[i] = c + end end diff --git a/src/core/constraint_template.jl b/src/core/constraint_template.jl index 6b7513cb..a2d94c23 100644 --- a/src/core/constraint_template.jl +++ b/src/core/constraint_template.jl @@ -28,154 +28,109 @@ end ### Nodal Constraints ### function constraint_flow_conservation(wm::AbstractWaterModel, i::Int; nw::Int=wm.cnw) # Collect various indices for edge-type components connected to node `i`. - check_valve_fr = _collect_comps_fr(wm, i, :check_valve; nw=nw) - check_valve_to = _collect_comps_to(wm, i, :check_valve; nw=nw) pipe_fr = _collect_comps_fr(wm, i, :pipe; nw=nw) pipe_to = _collect_comps_to(wm, i, :pipe; nw=nw) des_pipe_fr = _collect_comps_fr(wm, i, :des_pipe; nw=nw) des_pipe_to = _collect_comps_to(wm, i, :des_pipe; nw=nw) pump_fr = _collect_comps_fr(wm, i, :pump; nw=nw) pump_to = _collect_comps_to(wm, i, :pump; nw=nw) - pressure_reducing_valve_fr = _collect_comps_fr(wm, i, :pressure_reducing_valve; nw=nw) - pressure_reducing_valve_to = _collect_comps_to(wm, i, :pressure_reducing_valve; nw=nw) - shutoff_valve_fr = _collect_comps_fr(wm, i, :shutoff_valve; nw=nw) - shutoff_valve_to = _collect_comps_to(wm, i, :shutoff_valve; nw=nw) + regulator_fr = _collect_comps_fr(wm, i, :regulator; nw=nw) + regulator_to = _collect_comps_to(wm, i, :regulator; nw=nw) + short_pipe_fr = _collect_comps_fr(wm, i, :short_pipe; nw=nw) + short_pipe_to = _collect_comps_to(wm, i, :short_pipe; nw=nw) + valve_fr = _collect_comps_fr(wm, i, :valve; nw=nw) + valve_to = _collect_comps_to(wm, i, :valve; nw=nw) # Collect various indices for node-type components connected to node `i`. reservoirs = ref(wm, nw, :node_reservoir, i) # Reservoirs attached to node `i`. tanks = ref(wm, nw, :node_tank, i) # Tanks attached to node `i`. - junctions = ref(wm, nw, :node_junction, i) # Junctions attached to node `i`. + demands = ref(wm, nw, :node_demand, i) # Demands attached to node `i`. # Sum the constant demands required at node `i`. - nondispatchable_junctions = filter(j -> j in ids(wm, nw, :nondispatchable_junction), junctions) - fixed_demands = [ref(wm, nw, :nondispatchable_junction, j)["demand"] for j in nondispatchable_junctions] + nondispatchable_demands = filter(j -> j in ids(wm, nw, :nondispatchable_demand), demands) + fixed_demands = [ref(wm, nw, :nondispatchable_demand, j)["flow_rate"] for j in nondispatchable_demands] net_fixed_demand = length(fixed_demands) > 0 ? sum(fixed_demands) : 0.0 - # Get the indices of dispatchable junctions connected to node `i`. - dispatchable_junctions = filter(j -> j in ids(wm, nw, :dispatchable_junction), junctions) + # Get the indices of dispatchable demands connected to node `i`. + dispatchable_demands = filter(j -> j in ids(wm, nw, :dispatchable_demand), demands) # Initialize the flow conservation constraint dictionary entry. _initialize_con_dict(wm, :flow_conservation, nw=nw) # Add the flow conservation constraint. constraint_flow_conservation( - wm, nw, i, check_valve_fr, check_valve_to, pipe_fr, pipe_to, des_pipe_fr, - des_pipe_to, pump_fr, pump_to, pressure_reducing_valve_fr, - pressure_reducing_valve_to, shutoff_valve_fr, shutoff_valve_to, reservoirs, tanks, - dispatchable_junctions, net_fixed_demand) + wm, nw, i, pipe_fr, pipe_to, des_pipe_fr, des_pipe_to, pump_fr, pump_to, + regulator_fr, regulator_to, short_pipe_fr, short_pipe_to, valve_fr, valve_to, + reservoirs, tanks, dispatchable_demands, net_fixed_demand) end function constraint_node_directionality(wm::AbstractWaterModel, i::Int; nw::Int=wm.cnw) # Collect various indices for edge-type components connected to node `i`. - check_valve_fr = _collect_comps_fr(wm, i, :check_valve; nw=nw) - check_valve_to = _collect_comps_to(wm, i, :check_valve; nw=nw) pipe_fr = _collect_comps_fr(wm, i, :pipe; nw=nw) pipe_to = _collect_comps_to(wm, i, :pipe; nw=nw) des_pipe_fr = _collect_comps_fr(wm, i, :des_pipe; nw=nw) des_pipe_to = _collect_comps_to(wm, i, :des_pipe; nw=nw) pump_fr = _collect_comps_fr(wm, i, :pump; nw=nw) pump_to = _collect_comps_to(wm, i, :pump; nw=nw) - pressure_reducing_valve_fr = _collect_comps_fr(wm, i, :pressure_reducing_valve; nw=nw) - pressure_reducing_valve_to = _collect_comps_to(wm, i, :pressure_reducing_valve; nw=nw) - shutoff_valve_fr = _collect_comps_fr(wm, i, :shutoff_valve; nw=nw) - shutoff_valve_to = _collect_comps_to(wm, i, :shutoff_valve; nw=nw) + regulator_fr = _collect_comps_fr(wm, i, :regulator; nw=nw) + regulator_to = _collect_comps_to(wm, i, :regulator; nw=nw) + short_pipe_fr = _collect_comps_fr(wm, i, :short_pipe; nw=nw) + short_pipe_to = _collect_comps_to(wm, i, :short_pipe; nw=nw) + valve_fr = _collect_comps_fr(wm, i, :valve; nw=nw) + valve_to = _collect_comps_to(wm, i, :valve; nw=nw) + + # Collect various indices for node-type components connected to node `i`. + reservoirs = ref(wm, nw, :node_reservoir, i) # Reservoirs attached to node `i`. + tanks = ref(wm, nw, :node_tank, i) # Tanks attached to node `i`. + demands = ref(wm, nw, :node_demand, i) # Demands attached to node `i`. + + # Sum the constant demands required at node `i`. + nondispatchable_demands = filter(j -> j in ids(wm, nw, :nondispatchable_demand), demands) + fixed_demands = [ref(wm, nw, :nondispatchable_demand, j)["flow_rate"] for j in nondispatchable_demands] + net_fixed_demand = length(fixed_demands) > 0 ? sum(fixed_demands) : 0.0 # Get the number of nodal components attached to node `i`. - junctions = ref(wm, nw, :node_junction) - tanks = ref(wm, nw, :node_tank) - reservoirs = ref(wm, nw, :node_reservoir) - num_components = length(junctions) + length(tanks) + length(reservoirs) + num_components = length(demands) + length(tanks) + length(reservoirs) # Get the in degree of node `i`. - in_length = length(check_valve_to) + length(pipe_to) + length(des_pipe_to) + - length(pump_to) + length(pressure_reducing_valve_to) + length(shutoff_valve_to) + in_length = length(pipe_to) + length(des_pipe_to) + length(pump_to) + + length(regulator_to) + length(short_pipe_to) + length(valve_to) # Get the out degree of node `i`. - out_length = length(check_valve_fr) + length(pipe_fr) + length(des_pipe_fr) + - length(pump_fr) + length(pressure_reducing_valve_fr) + length(shutoff_valve_fr) + out_length = length(pipe_fr) + length(des_pipe_fr) + length(pump_fr) + + length(regulator_fr) + length(short_pipe_fr) + length(valve_fr) + + # Initialize the directionality constraint dictionary entry. + _initialize_con_dict(wm, :node_directionality, nw=nw) # Check if node directionality constraints should be added. if num_components == 0 && in_length + out_length == 2 - # Initialize the node directionality constraint dictionary entry. - _initialize_con_dict(wm, :node_directionality, nw=nw) - - # Add the node directionality constraint. - constraint_node_directionality( - wm, nw, i, check_valve_fr, check_valve_to, pipe_fr, pipe_to, des_pipe_fr, - des_pipe_to, pump_fr, pump_to, pressure_reducing_valve_fr, - pressure_reducing_valve_to, shutoff_valve_fr, shutoff_valve_to) - end -end - - -### Junction Constraints ### -function constraint_sink_directionality(wm::AbstractWaterModel, i::Int; nw::Int=wm.cnw) - # Collect various indices for edge-type components connected to node `i`. - check_valve_fr = _collect_comps_fr(wm, i, :check_valve; nw=nw) - check_valve_to = _collect_comps_to(wm, i, :check_valve; nw=nw) - pipe_fr = _collect_comps_fr(wm, i, :pipe; nw=nw) - pipe_to = _collect_comps_to(wm, i, :pipe; nw=nw) - des_pipe_fr = _collect_comps_fr(wm, i, :des_pipe; nw=nw) - des_pipe_to = _collect_comps_to(wm, i, :des_pipe; nw=nw) - pump_fr = _collect_comps_fr(wm, i, :pump; nw=nw) - pump_to = _collect_comps_to(wm, i, :pump; nw=nw) - pressure_reducing_valve_fr = _collect_comps_fr(wm, i, :pressure_reducing_valve; nw=nw) - pressure_reducing_valve_to = _collect_comps_to(wm, i, :pressure_reducing_valve; nw=nw) - shutoff_valve_fr = _collect_comps_fr(wm, i, :shutoff_valve; nw=nw) - shutoff_valve_to = _collect_comps_to(wm, i, :shutoff_valve; nw=nw) - - # Initialize the sink flow constraint dictionary entry. - _initialize_con_dict(wm, :sink_directionality, nw=nw) - - # Add the sink flow directionality constraint. - constraint_sink_directionality( - wm, nw, i, check_valve_fr, check_valve_to, pipe_fr, pipe_to, des_pipe_fr, - des_pipe_to, pump_fr, pump_to, pressure_reducing_valve_fr, - pressure_reducing_valve_to, shutoff_valve_fr, shutoff_valve_to) -end - - -### Reservoir Constraints ### -function constraint_reservoir_head(wm::AbstractWaterModel, i::Int; nw::Int=wm.cnw) - # Only fix the reservoir head if the reservoir is nondispatchable. - if !ref(wm, nw, :reservoir, i)["dispatchable"] - node_index = ref(wm, nw, :reservoir, i)["node"] - head = ref(wm, nw, :node, node_index)["head"] - _initialize_con_dict(wm, :reservoir_head, nw=nw) - constraint_reservoir_head(wm, nw, node_index, head) + # Add the intermediate node directionality constraint. + constraint_intermediate_directionality( + wm, nw, i, pipe_fr, pipe_to, des_pipe_fr, des_pipe_to, pump_fr, pump_to, + regulator_fr, regulator_to, short_pipe_fr, short_pipe_to, valve_fr, valve_to) + elseif length(reservoirs) > 0 && num_components == length(reservoirs) + # Add the source node directionality constraint. + constraint_source_directionality( + wm, nw, i, pipe_fr, pipe_to, des_pipe_fr, des_pipe_to, pump_fr, pump_to, + regulator_fr, regulator_to, short_pipe_fr, short_pipe_to, valve_fr, valve_to) + elseif num_components == length(fixed_demands) && net_fixed_demand < 0.0 + # Add the source node directionality constraint. + constraint_source_directionality( + wm, nw, i, pipe_fr, pipe_to, des_pipe_fr, des_pipe_to, pump_fr, pump_to, + regulator_fr, regulator_to, short_pipe_fr, short_pipe_to, valve_fr, valve_to) + elseif num_components == length(fixed_demands) && net_fixed_demand > 0.0 + # Add the sink node directionality constraint. + constraint_sink_directionality( + wm, nw, i, pipe_fr, pipe_to, des_pipe_fr, des_pipe_to, pump_fr, pump_to, + regulator_fr, regulator_to, short_pipe_fr, short_pipe_to, valve_fr, valve_to) end end -function constraint_source_directionality(wm::AbstractWaterModel, i::Int; nw::Int=wm.cnw) - # Collect various indices for edge-type components connected to node `i`. - check_valve_fr = _collect_comps_fr(wm, i, :check_valve; nw=nw) - check_valve_to = _collect_comps_to(wm, i, :check_valve; nw=nw) - pipe_fr = _collect_comps_fr(wm, i, :pipe; nw=nw) - pipe_to = _collect_comps_to(wm, i, :pipe; nw=nw) - des_pipe_fr = _collect_comps_fr(wm, i, :des_pipe; nw=nw) - des_pipe_to = _collect_comps_to(wm, i, :des_pipe; nw=nw) - pump_fr = _collect_comps_fr(wm, i, :pump; nw=nw) - pump_to = _collect_comps_to(wm, i, :pump; nw=nw) - pressure_reducing_valve_fr = _collect_comps_fr(wm, i, :pressure_reducing_valve; nw=nw) - pressure_reducing_valve_to = _collect_comps_to(wm, i, :pressure_reducing_valve; nw=nw) - shutoff_valve_fr = _collect_comps_fr(wm, i, :shutoff_valve; nw=nw) - shutoff_valve_to = _collect_comps_to(wm, i, :shutoff_valve; nw=nw) - - # Initialize the source flow constraint dictionary entry. - _initialize_con_dict(wm, :source_directionality, nw=nw) - - # Add the source flow directionality constraint. - constraint_source_directionality( - wm, nw, i, check_valve_fr, check_valve_to, pipe_fr, pipe_to, des_pipe_fr, - des_pipe_to, pump_fr, pump_to, pressure_reducing_valve_fr, - pressure_reducing_valve_to, shutoff_valve_fr, shutoff_valve_to) -end - - ### Tank Constraints ### -function constraint_tank_state(wm::AbstractWaterModel, i::Int; nw::Int=wm.cnw) +function constraint_tank_volume(wm::AbstractWaterModel, i::Int; nw::Int=wm.cnw) if !(:time_step in keys(ref(wm, nw))) Memento.error(_LOGGER, "Tank states cannot be controlled without a time step.") end @@ -186,13 +141,13 @@ function constraint_tank_state(wm::AbstractWaterModel, i::Int; nw::Int=wm.cnw) initial_level = tank["init_level"] surface_area = 0.25 * pi * tank["diameter"]^2 V_initial = surface_area * initial_level - _initialize_con_dict(wm, :tank_state, nw=nw) - constraint_tank_state_initial(wm, nw, i, V_initial) + _initialize_con_dict(wm, :tank_volume, nw=nw) + constraint_tank_volume_fixed(wm, nw, i, V_initial) end end -function constraint_tank_state(wm::AbstractWaterModel, i::Int, nw_1::Int, nw_2::Int) +function constraint_tank_volume(wm::AbstractWaterModel, i::Int, nw_1::Int, nw_2::Int) if !(:time_step in keys(ref(wm, nw_1))) Memento.error(_LOGGER, "Tank states cannot be controlled without a time step.") end @@ -201,113 +156,153 @@ function constraint_tank_state(wm::AbstractWaterModel, i::Int, nw_1::Int, nw_2:: if !ref(wm, nw_2, :tank, i)["dispatchable"] # TODO: What happens if a tank exists in nw_1 but not in nw_2? The index # "i" is assumed to be present in both when this constraint is applied. - _initialize_con_dict(wm, :tank_state, nw=nw_2) - constraint_tank_state(wm, nw_1, nw_2, i, ref(wm, nw_1, :time_step)) + _initialize_con_dict(wm, :tank_volume, nw=nw_2) + constraint_tank_volume(wm, nw_1, nw_2, i, ref(wm, nw_1, :time_step)) end end ### Pipe Constraints ### +function constraint_pipe_flow(wm::AbstractWaterModel, a::Int; nw::Int=wm.cnw, kwargs...) + _initialize_con_dict(wm, :pipe_flow, nw=nw, is_array=true) + con(wm, nw, :pipe_flow)[a] = Array{JuMP.ConstraintRef}([]) + constraint_pipe_flow(wm, nw, a) +end + + +function constraint_pipe_head(wm::AbstractWaterModel, a::Int; nw::Int=wm.cnw, kwargs...) + _initialize_con_dict(wm, :pipe_head, nw=nw, is_array=true) + con(wm, nw, :pipe_head)[a] = Array{JuMP.ConstraintRef}([]) + node_fr, node_to = ref(wm, nw, :pipe, a)["node_fr"], ref(wm, nw, :pipe, a)["node_to"] + constraint_pipe_head(wm, nw, a, node_fr, node_to) +end + + function constraint_pipe_head_loss(wm::AbstractWaterModel, a::Int; nw::Int=wm.cnw, kwargs...) + _initialize_con_dict(wm, :pipe_head_loss, nw=nw, is_array=true) + con(wm, nw, :pipe_head_loss)[a] = Array{JuMP.ConstraintRef}([]) node_fr, node_to = ref(wm, nw, :pipe, a)["node_fr"], ref(wm, nw, :pipe, a)["node_to"] - alpha, L = ref(wm, nw, :alpha), ref(wm, nw, :pipe, a)["length"] + exponent, L = ref(wm, nw, :alpha), ref(wm, nw, :pipe, a)["length"] r = minimum(ref(wm, nw, :resistance, a)) - - # Add common constraints used to model pump behavior. - _initialize_con_dict(wm, :pipe, nw=nw, is_array=true) - con(wm, nw, :pipe)[a] = Array{JuMP.ConstraintRef}([]) - constraint_pipe_common(wm, nw, a, node_fr, node_to, alpha, L, r) - - # Add constraints related to modeling head loss along pipe. - _initialize_con_dict(wm, :head_loss, nw=nw, is_array=true) - con(wm, nw, :head_loss)[a] = Array{JuMP.ConstraintRef}([]) - constraint_pipe_head_loss(wm, nw, a, node_fr, node_to, alpha, L, r) - constraint_pipe_head_loss_ub(wm, nw, a, alpha, L, r) + constraint_pipe_head_loss(wm, nw, a, node_fr, node_to, exponent, L, r) end -function constraint_pipe_head_loss_des(wm::AbstractWaterModel, a::Int; nw::Int=wm.cnw, kwargs...) - alpha = ref(wm, nw, :alpha) - pipe = ref(wm, nw, :des_pipe, a) +### Design Pipe Constraints ### +function constraint_on_off_pipe_flow_des(wm::AbstractWaterModel, a::Int; nw::Int=wm.cnw, kwargs...) + _initialize_con_dict(wm, :on_off_pipe_flow_des, nw=nw, is_array=true) + con(wm, nw, :on_off_pipe_flow_des)[a] = Array{JuMP.ConstraintRef}([]) resistances = ref(wm, nw, :resistance, a) - i, j = pipe["node_fr"], pipe["node_to"] + constraint_on_off_pipe_flow_des(wm, nw, a, resistances) +end + + +function constraint_on_off_pipe_head_des(wm::AbstractWaterModel, a::Int; nw::Int=wm.cnw, kwargs...) + _initialize_con_dict(wm, :on_off_pipe_head_des, nw=nw, is_array=true) + con(wm, nw, :on_off_pipe_head_des)[a] = Array{JuMP.ConstraintRef}([]) + node_fr = ref(wm, nw, :des_pipe, a)["node_fr"] + node_to = ref(wm, nw, :des_pipe, a)["node_to"] + constraint_on_off_pipe_head_des(wm, nw, a, node_fr, node_to) +end - _initialize_con_dict(wm, :head_loss, nw=nw, is_array=true) - con(wm, nw, :head_loss)[a] = Array{JuMP.ConstraintRef}([]) - constraint_flow_direction_selection_des(wm, a; nw=nw, kwargs...) - constraint_pipe_head_loss_des(wm, nw, a, alpha, i, j, pipe["length"], resistances) - constraint_pipe_head_loss_ub_des(wm, a; nw=nw, kwargs...) - constraint_resistance_selection_des(wm, a; nw=nw, kwargs...) +function constraint_on_off_pipe_head_loss_des(wm::AbstractWaterModel, a::Int; nw::Int=wm.cnw, kwargs...) + _initialize_con_dict(wm, :on_off_pipe_head_loss_des, nw=nw, is_array=true) + con(wm, nw, :on_off_pipe_head_loss_des)[a] = Array{JuMP.ConstraintRef}([]) + node_fr = ref(wm, nw, :des_pipe, a)["node_fr"] + node_to = ref(wm, nw, :des_pipe, a)["node_to"] + exponent, L = ref(wm, nw, :alpha), ref(wm, nw, :des_pipe, a)["length"] + resist = ref(wm, nw, :resistance, a) # Dictionary of possibel resistances at `a`. + constraint_on_off_pipe_head_loss_des(wm, nw, a, exponent, node_fr, node_to, L, resist) end -function constraint_resistance_selection_des(wm::AbstractWaterModel, a::Int; nw::Int=wm.cnw, kwargs...) - pipe_resistances = ref(wm, nw, :resistance, a) - constraint_resistance_selection_des(wm, nw, a, pipe_resistances; kwargs...) +### Pump Constraints ### +function constraint_on_off_pump_flow(wm::AbstractWaterModel, a::Int; nw::Int=wm.cnw, kwargs...) + node_fr, node_to = ref(wm, nw, :pump, a)["node_fr"], ref(wm, nw, :pump, a)["node_to"] + _initialize_con_dict(wm, :on_off_pump_flow, nw=nw, is_array=true) + con(wm, nw, :on_off_pump_flow)[a] = Array{JuMP.ConstraintRef}([]) + q_min_active = get(ref(wm, nw, :pump, a), "q_min_active", _q_eps) + constraint_on_off_pump_flow(wm, nw, a, q_min_active) end -function constraint_flow_direction_selection_des(wm::AbstractWaterModel, a::Int; nw::Int=wm.cnw) - pipe_resistances = ref(wm, nw, :resistance, a) - constraint_flow_direction_selection_des(wm, nw, a, pipe_resistances) +function constraint_on_off_pump_head(wm::AbstractWaterModel, a::Int; nw::Int=wm.cnw, kwargs...) + node_fr, node_to = ref(wm, nw, :pump, a)["node_fr"], ref(wm, nw, :pump, a)["node_to"] + _initialize_con_dict(wm, :on_off_pump_head, nw=nw, is_array=true) + con(wm, nw, :on_off_pump_head)[a] = Array{JuMP.ConstraintRef}([]) + constraint_on_off_pump_head(wm, nw, a, node_fr, node_to) end -function constraint_pipe_head_loss_ub_des(wm::AbstractWaterModel, a::Int; nw::Int=wm.cnw) - alpha, L = ref(wm, nw, :alpha), ref(wm, nw, :des_pipe, a)["length"] - pipe_resistances = ref(wm, nw, :resistance, a) - constraint_pipe_head_loss_ub_des(wm, nw, a, alpha, L, pipe_resistances) +function constraint_on_off_pump_head_gain(wm::AbstractWaterModel, a::Int; nw::Int=wm.cnw, kwargs...) + node_fr, node_to = ref(wm, nw, :pump, a)["node_fr"], ref(wm, nw, :pump, a)["node_to"] + _initialize_con_dict(wm, :on_off_pump_head_gain, nw=nw, is_array=true) + con(wm, nw, :on_off_pump_head_gain)[a] = Array{JuMP.ConstraintRef}([]) + coeffs = _get_function_from_head_curve(ref(wm, nw, :pump, a)["head_curve"]) + q_min_active = get(ref(wm, nw, :pump, a), "q_min_active", _q_eps) + constraint_on_off_pump_head_gain(wm, nw, a, node_fr, node_to, coeffs, q_min_active) end -### Check Valve Constraints ### -function constraint_check_valve_head_loss(wm::AbstractWaterModel, a::Int; nw::Int=wm.cnw, kwargs...) - node_fr = ref(wm, nw, :check_valve, a)["node_fr"] - node_to = ref(wm, nw, :check_valve, a)["node_to"] - alpha, L = ref(wm, nw, :alpha), ref(wm, nw, :check_valve, a)["length"] - r = maximum(ref(wm, nw, :resistance, a)) +### Short Pipe Constraints ### +function constraint_short_pipe_flow(wm::AbstractWaterModel, a::Int; nw::Int=wm.cnw, kwargs...) + _initialize_con_dict(wm, :short_pipe_flow, nw=nw, is_array=true) + con(wm, nw, :short_pipe_flow)[a] = Array{JuMP.ConstraintRef}([]) + constraint_short_pipe_flow(wm, nw, a) +end - # Add all common check valve constraints. - _initialize_con_dict(wm, :check_valve, nw=nw, is_array=true) - con(wm, nw, :check_valve)[a] = Array{JuMP.ConstraintRef}([]) - constraint_check_valve_common(wm, nw, a, node_fr, node_to) - # Add constraints that describe head loss relationships. - _initialize_con_dict(wm, :head_loss, nw=nw, is_array=true) - alpha, L = [ref(wm, nw, :alpha), ref(wm, nw, :check_valve, a)["length"]] - r = minimum(ref(wm, nw, :resistance, a)) - con(wm, nw, :head_loss)[a] = Array{JuMP.ConstraintRef}([]) - constraint_check_valve_head_loss(wm, nw, a, node_fr, node_to, L, r) - constraint_head_loss_ub_cv(wm, nw, a, alpha, L, r) +function constraint_short_pipe_head(wm::AbstractWaterModel, a::Int; nw::Int=wm.cnw, kwargs...) + node_fr = ref(wm, nw, :short_pipe, a)["node_fr"] + node_to = ref(wm, nw, :short_pipe, a)["node_to"] + _initialize_con_dict(wm, :short_pipe_head, nw=nw, is_array=true) + con(wm, nw, :short_pipe_head)[a] = Array{JuMP.ConstraintRef}([]) + constraint_short_pipe_head(wm, nw, a, node_fr, node_to) end -### Shutoff Valve Constraints ### -function constraint_shutoff_valve_head_loss(wm::AbstractWaterModel, a::Int; nw::Int=wm.cnw, kwargs...) - # Add all common shutoff valve constraints. - node_fr = ref(wm, nw, :shutoff_valve, a)["node_fr"] - node_to = ref(wm, nw, :shutoff_valve, a)["node_to"] - _initialize_con_dict(wm, :sv, nw=nw, is_array=true) - con(wm, nw, :sv)[a] = Array{JuMP.ConstraintRef}([]) - constraint_sv_common(wm, nw, a, node_fr, node_to) +### Valve Constraints ### +function constraint_on_off_valve_flow(wm::AbstractWaterModel, a::Int; nw::Int=wm.cnw, kwargs...) + node_fr, node_to = ref(wm, nw, :valve, a)["node_fr"], ref(wm, nw, :valve, a)["node_to"] + _initialize_con_dict(wm, :on_off_valve_flow, nw=nw, is_array=true) + con(wm, nw, :on_off_valve_flow)[a] = Array{JuMP.ConstraintRef}([]) + constraint_on_off_valve_flow(wm, nw, a) +end - # Add constraints that describe head loss relationships. - _initialize_con_dict(wm, :head_loss, nw=nw, is_array=true) - alpha, L = ref(wm, nw, :alpha), ref(wm, nw, :shutoff_valve, a)["length"] - r = minimum(ref(wm, nw, :resistance, a)) - con(wm, nw, :head_loss)[a] = Array{JuMP.ConstraintRef}([]) - constraint_shutoff_valve_head_loss(wm, nw, a, node_fr, node_to, L, r) - constraint_shutoff_valve_head_loss_ub(wm, nw, a, alpha, L, r) + +function constraint_on_off_valve_head(wm::AbstractWaterModel, a::Int; nw::Int=wm.cnw, kwargs...) + node_fr, node_to = ref(wm, nw, :valve, a)["node_fr"], ref(wm, nw, :valve, a)["node_to"] + _initialize_con_dict(wm, :on_off_valve_head, nw=nw, is_array=true) + con(wm, nw, :on_off_valve_head)[a] = Array{JuMP.ConstraintRef}([]) + constraint_on_off_valve_head(wm, nw, a, node_fr, node_to) +end + + +### Regulator Constraints ### +function constraint_on_off_regulator_flow(wm::AbstractWaterModel, a::Int; nw::Int=wm.cnw, kwargs...) + _initialize_con_dict(wm, :on_off_regulator_flow, nw=nw, is_array=true) + con(wm, nw, :on_off_regulator_flow)[a] = Array{JuMP.ConstraintRef}([]) + q_min_active = get(ref(wm, nw, :regulator, a), "q_min_active", _q_eps) + constraint_on_off_regulator_flow(wm, nw, a, q_min_active) +end + + +function constraint_on_off_regulator_head(wm::AbstractWaterModel, a::Int; nw::Int=wm.cnw, kwargs...) + _initialize_con_dict(wm, :on_off_regulator_head, nw=nw, is_array=true) + con(wm, nw, :on_off_regulator_head)[a] = Array{JuMP.ConstraintRef}([]) + node_fr = ref(wm, nw, :regulator, a)["node_fr"] + node_to = ref(wm, nw, :regulator, a)["node_to"] + elevation = ref(wm, nw, :node, node_to)["elevation"] + head_setting = elevation + ref(wm, nw, :regulator, a)["setting"] + constraint_on_off_regulator_head(wm, nw, a, node_fr, node_to, head_setting) end ### Pressure Reducing Valve Constraints ### function constraint_prv_head_loss(wm::AbstractWaterModel, a::Int; nw::Int=wm.cnw, kwargs...) - node_fr = ref(wm, nw, :pressure_reducing_valve, a)["node_fr"] - node_to = ref(wm, nw, :pressure_reducing_valve, a)["node_to"] - elev = ref(wm, nw, :node, node_to)["elevation"] - h_prv = elev + ref(wm, nw, :pressure_reducing_valve, a)["setting"] + node_fr = ref(wm, nw, :regulator, a)["node_fr"] + node_to = ref(wm, nw, :regulator, a)["node_to"] # Add all common pressure reducing valve constraints. _initialize_con_dict(wm, :prv, nw=nw, is_array=true) diff --git a/src/core/data.jl b/src/core/data.jl index 03ee25a8..6c9055f0 100644 --- a/src/core/data.jl +++ b/src/core/data.jl @@ -5,17 +5,60 @@ end """ -Turns given single network data into multinetwork data with a `count` replicate -of the given network. Note that this function performs a deepcopy of the -network data. Significant multinetwork space savings can often be achieved by -building application specific methods of building multinetwork with minimal -data replication. +Turns given single network data into multinetwork data with a `count` replicate of the given +network. Note that this function performs a deepcopy of the network data. Significant +multinetwork space savings can often be achieved by building application specific methods of +building multinetwork with minimal data replication. """ function replicate(data::Dict{String,<:Any}, count::Int; global_keys::Set{String}=Set{String}()) return _IM.replicate(data, count, union(global_keys, _wm_global_keys)) end +function _calc_pump_flow_min(pump::Dict{String,<:Any}, node_fr::Dict{String,Any}, node_to::Dict{String,Any}) + return get(pump, "flow_min", 0.0) +end + + +function _calc_pump_flow_min_forward(pump::Dict{String,<:Any}, node_fr::Dict{String,Any}, node_to::Dict{String,Any}) + flow_min_forward = get(pump, "flow_min_forward", _q_eps) + return max(_calc_pump_flow_min(pump, node_fr, node_to), flow_min_forward) +end + + +function _calc_pump_flow_max_reverse(pump::Dict{String,<:Any}, node_fr::Dict{String,Any}, node_to::Dict{String,Any}) + flow_max_reverse = get(pump, "flow_max_reverse", 0.0) + return min(_calc_pump_flow_max(pump, node_fr, node_to), flow_max_reverse) +end + + +function _calc_pump_flow_max(pump::Dict{String,<:Any}, node_fr::Dict{String,Any}, node_to::Dict{String,Any}) + coeff = _get_function_from_head_curve(pump["head_curve"]) + q_max_1 = (-coeff[2] + sqrt(coeff[2]^2 - 4.0*coeff[1]*coeff[3])) * inv(2.0*coeff[1]) + q_max_2 = (-coeff[2] - sqrt(coeff[2]^2 - 4.0*coeff[1]*coeff[3])) * inv(2.0*coeff[1]) + return min(max(q_max_1, q_max_2), get(pump, "flow_max", Inf)) +end + + +function _calc_pump_flow_bounds_active(pump::Dict{String,<:Any}) + q_min, q_max = _calc_pump_flow_bounds(pump) + q_min = max(max(get(pump, "flow_min_forward", _q_eps), _q_eps), q_min) + return q_min, q_max +end + + +function _correct_flow_bounds!(data::Dict{String,<:Any}) + for (idx, pump) in get(data, "pump", Dict{String,Any}()) + node_fr_id, node_to_id = string(pump["node_fr"]), string(pump["node_to"]) + node_fr, node_to = data["node"][node_fr_id], data["node"][node_to_id] + pump["flow_min"] = _calc_pump_flow_min(pump, node_fr, node_to) + pump["flow_max"] = _calc_pump_flow_max(pump, node_fr, node_to) + pump["flow_min_forward"] = _calc_pump_flow_min_forward(pump, node_fr, node_to) + pump["flow_max_reverse"] = _calc_pump_flow_max_reverse(pump, node_fr, node_to) + end +end + + function calc_resistances_hw(pipes::Dict{<:Any, <:Any}) resistances = Dict([(a, Array{Float64, 1}()) for a in keys(pipes)]) @@ -206,7 +249,7 @@ end function set_start_reservoir!(data) for (i, reservoir) in data["reservoir"] - reservoir["qr_start"] = reservoir["qr"] + reservoir["q_reservoir_start"] = reservoir["q"] end end diff --git a/src/core/objective.jl b/src/core/objective.jl index 28106d02..c8c3167b 100644 --- a/src/core/objective.jl +++ b/src/core/objective.jl @@ -16,18 +16,18 @@ end """ objective_des(wm::AbstractWaterModel) -Sets the objective function for network design (des) problem specifications. -By default, the cost of selecting discrete network resistances is minimized. +Sets the objective function for network design (des) problem specifications. By default, the +cost of selecting the discrete pipe resistances over all design pipes is minimized. """ function objective_des(wm::AbstractWaterModel) - obj = JuMP.AffExpr(0.0) + objective = JuMP.AffExpr(0.0) - for n in nw_ids(wm) - for (a, pipe) in ref(wm, n, :des_pipe) - costs = pipe["length"] .* ref(wm, n, :resistance_cost, a) - JuMP.add_to_expression!(obj, sum(costs .* var(wm, n, :x_res, a))) + for (n, nw_ref) in nws(wm) + for (a, des_pipe) in ref(wm, n, :des_pipe) + costs = des_pipe["length"] .* ref(wm, n, :resistance_cost, a) + JuMP.add_to_expression!(objective, sum(costs .* var(wm, n, :z_des_pipe, a))) end end - JuMP.@objective(wm.model, _MOI.MIN_SENSE, obj) + JuMP.@objective(wm.model, _MOI.MIN_SENSE, objective) end diff --git a/src/core/ref.jl b/src/core/ref.jl index 3fed696c..495e41c1 100644 --- a/src/core/ref.jl +++ b/src/core/ref.jl @@ -33,12 +33,12 @@ end function calc_demand_bounds(wm::AbstractWaterModel, n::Int=wm.cnw) # Initialize the dictionaries for minimum and maximum heads. - demand_min = Dict((i, -Inf) for (i, junc) in ref(wm, n, :dispatchable_junction)) - demand_max = Dict((i, Inf) for (i, junc) in ref(wm, n, :dispatchable_junction)) + demand_min = Dict((i, -Inf) for (i, demand) in ref(wm, n, :dispatchable_demand)) + demand_max = Dict((i, Inf) for (i, demand) in ref(wm, n, :dispatchable_demand)) - for (i, junction) in ref(wm, n, :dispatchable_junction) - demand_min[i] = max(demand_min[i], junction["demand_min"]) - demand_max[i] = min(demand_max[i], junction["demand_max"]) + for (i, demand) in ref(wm, n, :dispatchable_demand) + demand_min[i] = max(demand_min[i], demand["demand_min"]) + demand_max[i] = min(demand_max[i], demand["demand_max"]) end return demand_min, demand_max @@ -66,12 +66,12 @@ function calc_head_bounds(wm::AbstractWaterModel, n::Int=wm.cnw) head_max = Dict((i, max_head) for (i, node) in ref(wm, n, :node)) for (i, node) in ref(wm, n, :node) - num_junctions = length(ref(wm, n, :node_junction, i)) + num_demands = length(ref(wm, n, :node_demand, i)) num_reservoirs = length(ref(wm, n, :node_reservoir, i)) num_tanks = length(ref(wm, n, :node_tank, i)) # If the node has zero demand, pressures can be negative. - if num_junctions + num_reservoirs + num_tanks == 0 + if num_demands + num_reservoirs + num_tanks == 0 head_min[i] -= 100.0 end end @@ -96,9 +96,9 @@ function calc_head_bounds(wm::AbstractWaterModel, n::Int=wm.cnw) head_max[node_id] = elevation + tank["max_level"] end - for (a, pressure_reducing_valve) in ref(wm, n, :pressure_reducing_valve) - p_setting = ref(wm, n, :pressure_reducing_valve, a)["setting"] - node_to = pressure_reducing_valve["node_to"] + for (a, regulator) in ref(wm, n, :regulator) + p_setting = ref(wm, n, :regulator, a)["setting"] + node_to = regulator["node_to"] h_setting = ref(wm, n, :node, node_to)["elevation"] + p_setting head_min[node_to] = min(head_min[node_to], h_setting) end @@ -117,10 +117,10 @@ function calc_flow_bounds(wm::AbstractWaterModel, n::Int=wm.cnw) h_lb, h_ub = calc_head_bounds(wm, n) alpha = ref(wm, n, :alpha) - nondispatchable_junctions = values(ref(wm, n, :nondispatchable_junction)) - sum_demand = length(nondispatchable_junctions) > 0 ? sum(junc["demand"] for junc in nondispatchable_junctions) : 0.0 - dispatchable_junctions = values(ref(wm, n, :dispatchable_junction)) - sum_demand += length(dispatchable_junctions) > 0 ? sum(junc["demand_max"] for junc in dispatchable_junctions) : 0.0 + nondispatchable_demands = values(ref(wm, n, :nondispatchable_demand)) + sum_demand = length(nondispatchable_demands) > 0 ? sum(demand["flow_rate"] for demand in nondispatchable_demands) : 0.0 + dispatchable_demands = values(ref(wm, n, :dispatchable_demand)) + sum_demand += length(dispatchable_demands) > 0 ? sum(demand["demand_max"] for demand in dispatchable_demands) : 0.0 if :time_step in keys(ref(wm, n)) time_step = ref(wm, n, :time_step) @@ -130,7 +130,7 @@ function calc_flow_bounds(wm::AbstractWaterModel, n::Int=wm.cnw) lb, ub = Dict{String,Any}(), Dict{String,Any}() - for name in ["check_valve", "pipe", "des_pipe", "shutoff_valve"] + for name in ["pipe", "des_pipe"] lb[name], ub[name] = Dict{Int,Array{Float64}}(), Dict{Int,Array{Float64}}() for (a, comp) in ref(wm, n, Symbol(name)) @@ -167,14 +167,38 @@ function calc_flow_bounds(wm::AbstractWaterModel, n::Int=wm.cnw) end end - lb["pressure_reducing_valve"] = Dict{Int,Array{Float64}}() - ub["pressure_reducing_valve"] = Dict{Int,Array{Float64}}() + lb["regulator"] = Dict{Int,Array{Float64}}() + ub["regulator"] = Dict{Int,Array{Float64}}() - for (a, prv) in ref(wm, n, :pressure_reducing_valve) - name = "pressure_reducing_valve" + for (a, regulator) in ref(wm, n, :regulator) + name = "regulator" lb[name][a], ub[name][a] = [0.0], [sum_demand] - haskey(prv, "q_min") && (lb[name][a][1] = max(lb[name][a][1], prv["q_min"])) - haskey(prv, "q_max") && (ub[name][a][1] = min(ub[name][a][1], prv["q_max"])) + haskey(regulator, "q_min") && (lb[name][a][1] = max(lb[name][a][1], regulator["q_min"])) + haskey(regulator, "q_max") && (ub[name][a][1] = min(ub[name][a][1], regulator["q_max"])) + end + + lb["short_pipe"] = Dict{Int,Float64}() + ub["short_pipe"] = Dict{Int,Float64}() + + for (a, short_pipe) in ref(wm, n, :short_pipe) + lb["short_pipe"][a], ub["short_pipe"][a] = -sum_demand, sum_demand + haskey(short_pipe, "q_min") && (lb["short_pipe"][a] = max(lb["short_pipe"][a], short_pipe["q_min"])) + haskey(short_pipe, "q_max") && (ub["short_pipe"][a] = min(ub["short_pipe"][a], short_pipe["q_max"])) + + if short_pipe["flow_direction"] == POSITIVE + lb["short_pipe"][a] = max(lb["short_pipe"][a], 0.0) + elseif short_pipe["flow_direction"] == NEGATIVE + ub["short_pipe"][a] = min(ub["short_pipe"][a], 0.0) + end + end + + lb["valve"] = Dict{Int,Array{Float64}}() + ub["valve"] = Dict{Int,Array{Float64}}() + + for (a, valve) in ref(wm, n, :valve) + lb["valve"][a], ub["valve"][a] = [-sum_demand], [sum_demand] + haskey(valve, "q_min") && (lb["valve"][a][1] = max(lb["valve"][a][1], valve["q_min"])) + haskey(valve, "q_max") && (ub["valve"][a][1] = min(ub["valve"][a][1], valve["q_max"])) end lb["pump"], ub["pump"] = Dict{Int,Array{Float64}}(), Dict{Int,Array{Float64}}() diff --git a/src/core/types.jl b/src/core/types.jl index 8898422d..c3502151 100644 --- a/src/core/types.jl +++ b/src/core/types.jl @@ -8,13 +8,15 @@ const _q_eps = 6.31465679e-6 "Models derived from AbstractDirectedFlowModel" mutable struct CRDWaterModel <: AbstractWaterModel @wm_fields end mutable struct LRDWaterModel <: AbstractWaterModel @wm_fields end +mutable struct QRDWaterModel <: AbstractWaterModel @wm_fields end +mutable struct CQRDWaterModel <: AbstractWaterModel @wm_fields end "Models derived from AbstractUndirectedFlowModel" mutable struct NCWaterModel <: AbstractWaterModel @wm_fields end mutable struct LAWaterModel <: AbstractWaterModel @wm_fields end "Models that use direction variables." -AbstractDirectedModel = Union{CRDWaterModel,LRDWaterModel} +AbstractDirectedModel = Union{CRDWaterModel,LRDWaterModel,QRDWaterModel,CQRDWaterModel} "Models that don't use direction variables." AbstractUndirectedModel = Union{NCWaterModel,LAWaterModel} diff --git a/src/core/variable.jl b/src/core/variable.jl index 47737c2d..01044b72 100644 --- a/src/core/variable.jl +++ b/src/core/variable.jl @@ -2,16 +2,19 @@ # This file defines commonly-used variables for water systems models. ######################################################################## + "Sets the start value for a given variable." function comp_start_value(comp::Dict{String,<:Any}, key::String, default::Float64=0.0) return get(comp, key, default) end + "Sets the start value for a given variable." function comp_start_value(comp::Dict{String,<:Any}, key_1::String, key_2::Int64, default=0.0) return key_1 in keys(comp) ? get(get(comp, key_1, default), key_2, default) : default end + "Given a variable that is indexed by component IDs, builds the standard solution structure." function sol_component_value(wm::AbstractWaterModel, n::Int, comp_name::Symbol, field_name::Symbol, comp_ids, variables) for i in comp_ids @@ -20,153 +23,139 @@ function sol_component_value(wm::AbstractWaterModel, n::Int, comp_name::Symbol, end end + ### Variables related to nodal components. ### "Creates bounded (by default) or unbounded total hydraulic head (or head) variables for all nodes in the network, i.e., `h[i]` for `i` in `node`." function variable_head(wm::AbstractWaterModel; nw::Int=wm.cnw, bounded::Bool=true, report::Bool=true) + # Initialize variables for total hydraulic head. h = var(wm, nw)[:h] = JuMP.@variable(wm.model, [i in ids(wm, nw, :node)], base_name="$(nw)_h", start=comp_start_value(ref(wm, nw, :node, i), "h_start")) if bounded + # Compute the lower and upper head bound dictionaries. h_lb, h_ub = calc_head_bounds(wm, nw) for (i, node) in ref(wm, nw, :node) + # Set the lower and upper bounds for each head. JuMP.set_lower_bound(h[i], h_lb[i]) JuMP.set_upper_bound(h[i], h_ub[i]) end end + # Initialize an entry to the solution component dictionary for heads. report && sol_component_value(wm, nw, :node, :h, ids(wm, nw, :node), h) - # Create expressions that calculate pressures. + # Create expressions that calculate pressures from head. p = var(wm, nw)[:p] = JuMP.@expression(wm.model, - [i in ids(wm, nw, :node)], - var(wm, nw, :h, i) - ref(wm, nw, :node, i)["elevation"]) + [i in ids(wm, nw, :node)], var(wm, nw, :h, i) - ref(wm, nw, :node, i)["elevation"]) + # Initialize an entry to the solution component dictionary for pressures. report && sol_component_value(wm, nw, :node, :p, ids(wm, nw, :node), p) + + # Create an expression that maps total hydraulic head to volume for tanks. + V = var(wm, nw)[:V] = JuMP.@expression(wm.model, [i in ids(wm, nw, :tank)], + (var(wm, nw, :h, ref(wm, nw, :tank, i)["node"]) - + ref(wm, nw, :node, ref(wm, nw, :tank, i)["node"])["elevation"]) + * 0.25 * pi * ref(wm, nw, :tank, i)["diameter"]^2) + + # Initialize an entry to the solution component dictionary for volumes. + report && sol_component_value(wm, nw, :tank, :V, ids(wm, nw, :tank), V) end -"Creates head gain variables corresponding to all pumps in the network, i.e., -`g[a]` for `a` in `pump`. These denote head gains between nodes i and j." -function variable_head_gain(wm::AbstractWaterModel; nw::Int=wm.cnw, report::Bool=true) - g = var(wm, nw)[:g] = JuMP.@variable(wm.model, [a in ids(wm, nw, :pump)], - base_name="$(nw)_g", lower_bound=0.0, # Pump flow is nonnegative. - start=comp_start_value(ref(wm, nw, :pump, a), "g_start")) +function variable_pump_head_gain(wm::AbstractWaterModel; nw::Int=wm.cnw, bounded::Bool=true, report::Bool=true) + # Initialize variables for total hydraulic head gain from a pump. + g = var(wm, nw)[:g_pump] = JuMP.@variable(wm.model, [a in ids(wm, nw, :pump)], + base_name="$(nw)_g_pump", lower_bound=0.0, # Pump gain is nonnegative. + start=comp_start_value(ref(wm, nw, :pump, a), "g_pump_start", 1.0e-6)) + + # Initialize an entry to the solution component dictionary for head gains. report && sol_component_value(wm, nw, :pump, :g, ids(wm, nw, :pump), g) end -"Creates outgoing flow variables for all reservoirs in the network, i.e., -`qr[i]` for `i` in `reservoir`. Note that these variables are always -nonnegative, as there is never incoming flow to a reservoir." -function variable_reservoir(wm::AbstractWaterModel; nw::Int=wm.cnw, report::Bool=true) - qr = var(wm, nw)[:qr] = JuMP.@variable(wm.model, - [i in ids(wm, nw, :reservoir)], lower_bound=0.0, base_name="$(nw)_qr", - start=comp_start_value(ref(wm, nw, :reservoir, i), "qr_start")) - report && sol_component_value(wm, nw, :reservoir, :qr, ids(wm, nw, :reservoir), qr) +"Instantiates outgoing flow variables for all reservoirs in the network, i.e., +`q_reservoir[i]` for `i` in `reservoir`. Note that these variables are always nonnegative, +since for each reservoir, there will never be incoming flow." +function variable_reservoir_flow(wm::AbstractWaterModel; nw::Int=wm.cnw, report::Bool=true) + q_reservoir = var(wm, nw)[:q_reservoir] = JuMP.@variable(wm.model, + [i in ids(wm, nw, :reservoir)], lower_bound=0.0, base_name="$(nw)_q_reservoir", + start=comp_start_value(ref(wm, nw, :reservoir, i), "q_reservoir_start")) + + report && sol_component_value(wm, nw, :reservoir, :q, ids(wm, nw, :reservoir), q_reservoir) end -"Creates demand variables for all dispatchable junctions in the network, i.e., `demand[i]` -for `i` in `dispatchable_junction`." -function variable_dispatchable_junction(wm::AbstractWaterModel; nw::Int=wm.cnw, bounded::Bool=true, report::Bool=true) - demand = var(wm, nw)[:demand] = JuMP.@variable(wm.model, - [i in ids(wm, nw, :dispatchable_junction)], base_name="$(nw)_demand", - start=comp_start_value(ref(wm, nw, :dispatchable_junction, i), "demand_start")) + +"Instantiates demand flow variables for all dispatchable demands in the network, i.e., +`demand[i]` for `i` in `dispatchable_demand`." +function variable_demand_flow(wm::AbstractWaterModel; nw::Int=wm.cnw, bounded::Bool=true, report::Bool=true) + q_demand = var(wm, nw)[:q_demand] = JuMP.@variable(wm.model, + [i in ids(wm, nw, :dispatchable_demand)], base_name="$(nw)_q_demand", + start=comp_start_value(ref(wm, nw, :dispatchable_demand, i), "q_demand_start")) if bounded - demand_lb, demand_ub = calc_demand_bounds(wm, nw) + q_demand_lb, q_demand_ub = calc_demand_bounds(wm, nw) - for (i, junction) in ref(wm, nw, :dispatchable_junction) - JuMP.set_lower_bound(demand[i], demand_lb[i]) - JuMP.set_upper_bound(demand[i], demand_ub[i]) + for (i, demand) in ref(wm, nw, :dispatchable_demand) + JuMP.set_lower_bound(q_demand[i], q_demand_lb[i]) + JuMP.set_upper_bound(q_demand[i], q_demand_ub[i]) end end - report && sol_component_value(wm, nw, :junction, :demand, ids(wm, nw, :dispatchable_junction), demand) + report && sol_component_value(wm, nw, :demand, :q, ids(wm, nw, :dispatchable_demand), q_demand) end -"Creates outgoing flow variables for all tanks in the network, i.e., `qt[i]` -for `i` in `tank`. Note that, unlike reservoirs, tanks can have inflow." -function variable_tank(wm::AbstractWaterModel; nw::Int=wm.cnw, report::Bool=true) - qt = var(wm, nw)[:qt] = JuMP.@variable(wm.model, - [i in ids(wm, nw, :tank)], base_name="$(nw)_qt", - start=comp_start_value(ref(wm, nw, :tank, i), "qt_start")) - - report && sol_component_value(wm, nw, :tank, :qt, ids(wm, nw, :tank), qt) - # Create an expression that maps total hydraulic head to volume. - V = var(wm, nw)[:V] = JuMP.@expression(wm.model, - [i in ids(wm, nw, :tank)], - (var(wm, nw, :h, ref(wm, nw, :tank, i)["node"]) - - ref(wm, nw, :node, ref(wm, nw, :tank, i)["node"])["elevation"]) - * 0.25 * pi * ref(wm, nw, :tank, i)["diameter"]^2) +"Creates outgoing flow variables for all tanks in the network, i.e., `q_tank[i]` +for `i` in `tank`. Note that, unlike reservoirs, tanks can have inflow." +function variable_tank_flow(wm::AbstractWaterModel; nw::Int=wm.cnw, report::Bool=true) + q_tank = var(wm, nw)[:q_tank] = JuMP.@variable(wm.model, + [i in ids(wm, nw, :tank)], base_name="$(nw)_q_tank", + start=comp_start_value(ref(wm, nw, :tank, i), "q_tank_start")) - report && sol_component_value(wm, nw, :tank, :V, ids(wm, nw, :tank), V) + report && sol_component_value(wm, nw, :tank, :q, ids(wm, nw, :tank), q_tank) end ### Link variables. ### -"Creates binary variables for all check valves in the network, i.e., -`z_check_valve[a]` for `a` in `check_valve`, where one denotes that the check -valve is open and zero denotes that the check valve is closed." -function variable_check_valve_indicator(wm::AbstractWaterModel; nw::Int=wm.cnw, relax::Bool=false, report::Bool=true) +"Creates binary variables for valves in the network, i.e., `z_valve[a]` for `a` in `valve`, +where one denotes that the valve is open and zero denotes that the valve is closed." +function variable_valve_indicator(wm::AbstractWaterModel; nw::Int=wm.cnw, relax::Bool=false, report::Bool=true) if !relax - z_check_valve = var(wm, nw)[:z_check_valve] = JuMP.@variable(wm.model, - [a in ids(wm, nw, :check_valve)], base_name = "$(nw)_z_check_valve", - binary = true, - start = comp_start_value(ref(wm, nw, :check_valve, a), "z_check_valve_start")) + z_valve = var(wm, nw)[:z_valve] = JuMP.@variable(wm.model, + [a in ids(wm, nw, :valve)], base_name = "$(nw)_z_valve", binary = true, + start = comp_start_value(ref(wm, nw, :valve, a), "z_valve_start")) else - z_check_valve = var(wm, nw)[:z_check_valve] = JuMP.@variable(wm.model, - [a in ids(wm, nw, :check_valve)], base_name = "z_check_valve[$(nw)]", + z_valve = var(wm, nw)[:z_valve] = JuMP.@variable(wm.model, + [a in ids(wm, nw, :valve)], base_name = "z_valve[$(nw)]", lower_bound = 0.0, upper_bound = 1.0, - start = comp_start_value(ref(wm, nw, :check_valve, a), "z_check_valve_start")) + start = comp_start_value(ref(wm, nw, :valve, a), "z_valve_start")) end - report && _IM.sol_component_value(wm, nw, :check_valve, :status, ids(wm, nw, :check_valve), z_check_valve) + report && _IM.sol_component_value(wm, nw, :valve, :status, ids(wm, nw, :valve), z_valve) end -"Creates binary variables for all shutoff valves in the network, i.e., -`z_shutoff_valve[a]` for `a` in `sv`, where one denotes that the shutoff valve -is open and zero denotes that the shutoff valve is closed." -function variable_shutoff_valve_indicator(wm::AbstractWaterModel; nw::Int=wm.cnw, relax::Bool=false, report::Bool=true) - if !relax - z_shutoff_valve = var(wm, nw)[:z_shutoff_valve] = JuMP.@variable(wm.model, - [a in ids(wm, nw, :shutoff_valve)], base_name = "$(nw)_z_shutoff_valve", - binary = true, - start = comp_start_value(ref(wm, nw, :shutoff_valve, a), "z_shutoff_valve_start")) - else - z_shutoff_valve = var(wm, nw)[:z_shutoff_valve] = JuMP.@variable(wm.model, - [a in ids(wm, nw, :shutoff_valve)], base_name = "$(nw)_z_shutoff_valve", - lower_bound = 0.0, upper_bound = 1.0, - start = comp_start_value(ref(wm, nw, :shutoff_valve, a), "z_shutoff_valve_start")) - end - - report && sol_component_value(wm, nw, :shutoff_valve, :status, ids(wm, nw, :shutoff_valve), z_shutoff_valve) -end -"Creates binary variables for all PRVs in the network, i.e., -`z_pressure_reducing_valve[a]` for `a` in `pressure_reducing_valve`, where one -denotes that the pressure reducing is currently open and zero otherwise." -function variable_pressure_reducing_valve_indicator(wm::AbstractWaterModel; nw::Int=wm.cnw, relax::Bool=false, report::Bool=true) +"Creates binary variables for all regulators in the network, i.e., `z_regulator[a]` for `a` in +`regulator`, where one denotes that the pressure reducing is currently open and zero otherwise." +function variable_regulator_indicator(wm::AbstractWaterModel; nw::Int=wm.cnw, relax::Bool=false, report::Bool=true) if !relax - z_pressure_reducing_valve = var(wm, nw)[:z_pressure_reducing_valve] = JuMP.@variable(wm.model, - [a in ids(wm, nw, :pressure_reducing_valve)], base_name = - "$(nw)_z_pressure_reducing_valve", binary = true, start = - comp_start_value(ref(wm, nw, :pressure_reducing_valve, a), - "z_pressure_reducing_valve_start")) + z_regulator = var(wm, nw)[:z_regulator] = JuMP.@variable(wm.model, + [a in ids(wm, nw, :regulator)], base_name = "$(nw)_z_regulator", binary = true, + start = comp_start_value(ref(wm, nw, :regulator, a), "z_regulator_start")) else - z_pressure_reducing_valve = var(wm, nw)[:z_pressure_reducing_valve] = JuMP.@variable(wm.model, - [a in ids(wm, nw, :pressure_reducing_valve)], base_name = "$(nw)_z_pressure_reducing_valve", + z_regulator = var(wm, nw)[:z_regulator] = JuMP.@variable(wm.model, + [a in ids(wm, nw, :regulator)], base_name = "$(nw)_z_regulator", lower_bound = 0.0, upper_bound = 1.0, - start = comp_start_value(ref(wm, nw, :pressure_reducing_valve, a), "z_pressure_reducing_valve_start")) + start = comp_start_value(ref(wm, nw, :regulator, a), "z_regulator_start")) end - report && sol_component_value(wm, nw, :pressure_reducing_valve, :z_pressure_reducing_valve, - ids(wm, nw, :pressure_reducing_valve), z_pressure_reducing_valve) + report && sol_component_value(wm, nw, :regulator, :status, ids(wm, nw, :regulator), z_regulator) end + "Creates binary variables for all pumps in the network, i.e., `z_pump[a]` for `a` in `pump`, where one denotes that the pump is currently operating (i.e., on), and zero indicates that the pump is not operating (i.e., off)." @@ -174,8 +163,7 @@ function variable_pump_indicator(wm::AbstractWaterModel; nw::Int=wm.cnw, relax:: if !relax z_pump = var(wm, nw)[:z_pump] = JuMP.@variable(wm.model, [a in ids(wm, nw, :pump)], base_name = "$(nw)_z_pump", - binary = true, - start = comp_start_value(ref(wm, nw, :pump, a), "z_pump_start")) + binary = true, start = comp_start_value(ref(wm, nw, :pump, a), "z_pump_start")) else z_pump = var(wm, nw)[:z_pump] = JuMP.@variable(wm.model, [a in ids(wm, nw, :pump)], base_name = "$(nw)_z_pump", @@ -186,25 +174,26 @@ function variable_pump_indicator(wm::AbstractWaterModel; nw::Int=wm.cnw, relax:: report && sol_component_value(wm, nw, :pump, :status, ids(wm, nw, :pump), z_pump) end -function variable_pressure_reducing_valve_operation(wm::AbstractWaterModel; nw::Int=wm.cnw, report::Bool=true) - variable_pressure_reducing_valve_indicator(wm, nw=nw, report=report) -end - -function variable_pump_operation(wm::AbstractWaterModel; nw::Int=wm.cnw, report::Bool=true) -end -"Creates binary variables for all network design or design resistances in -the network, i.e., `x_res[a]` for `a` in `pipe`, for `r` in `resistance[a]`, -where one denotes that the given resistance is active in the design." -function variable_resistance(wm::AbstractWaterModel; nw::Int=wm.cnw, report::Bool=true) - x_res = var(wm, nw)[:x_res] = Dict{Int, Array{JuMP.VariableRef}}() +"Creates binary variables for all network design or design resistances in the network, i.e., +`z_des_pipe[a]` for `a` in `pipe`, for `r` in `resistance[a]`, where one denotes that the +given resistance is active in the design." +function variable_pipe_des_indicator(wm::AbstractWaterModel; nw::Int=wm.cnw, relax::Bool=false, report::Bool=true) + z = var(wm, nw)[:z_des_pipe] = Dict{Int, Array{JuMP.VariableRef}}() for a in ids(wm, nw, :des_pipe) n_r = length(ref(wm, nw, :resistance, a)) # Number of resistances. - var(wm, nw, :x_res)[a] = JuMP.@variable(wm.model, [r in 1:n_r], - binary=true, base_name="$(nw)_x_res[$(a)]", - start=comp_start_value(ref(wm, nw, :des_pipe, a), "x_res_start", r)) + + if !relax + var(wm, nw, :z_des_pipe)[a] = JuMP.@variable(wm.model, [r in 1:n_r], + binary=true, base_name="$(nw)_z_des_pipe[$(a)]", + start=comp_start_value(ref(wm, nw, :des_pipe, a), "z_des_pipe_start", r)) + else + var(wm, nw, :z_des_pipe)[a] = JuMP.@variable(wm.model, [r in 1:n_r], + lower_bound = 0.0, upper_bound = 0.0, base_name="$(nw)_z_des_pipe[$(a)]", + start=comp_start_value(ref(wm, nw, :des_pipe, a), "z_des_pipe_start", r)) + end end - report && sol_component_value(wm, nw, :des_pipe, :x_res, ids(wm, nw, :des_pipe), x_res) + report && sol_component_value(wm, nw, :des_pipe, :built, ids(wm, nw, :des_pipe), z) end diff --git a/src/form/cqrd.jl b/src/form/cqrd.jl new file mode 100644 index 00000000..1a7c92b6 --- /dev/null +++ b/src/form/cqrd.jl @@ -0,0 +1,124 @@ +# Define common CQRD (convex quadratic relaxation- and direction-based) implementations of +# water distribution network constraints, which use directed flow variables. + + +function constraint_pipe_head_loss(wm::CQRDWaterModel, n::Int, a::Int, node_fr::Int, node_to::Int, exponent::Float64, L::Float64, r::Float64) + # Get the number of breakpoints for the pipe. + num_breakpoints = get(wm.ext, :pipe_breakpoints, 1) + + # 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) + + # Loop over breakpoints strictly between the lower and upper variable bounds. + for pt in range(0.0, stop = JuMP.upper_bound(qp), length = num_breakpoints+2)[2:end-1] + # Add a linear outer approximation of the convex relaxation at `pt`. + lhs = _get_head_loss_oa_binary(qp, y, pt, exponent) + c = JuMP.@constraint(wm.model, r * lhs <= inv(L) * dhp) + append!(con(wm, n, :pipe_head_loss)[a], [c]) + end + + # Add linear upper bounds on the above outer approximations. + rhs = r * JuMP.upper_bound(qp)^(exponent - 1.0) * qp + c = JuMP.@constraint(wm.model, inv(L) * dhp <= rhs) + + # Append the :pipe_head_loss constraint array. + append!(con(wm, n, :pipe_head_loss)[a], [c]) + + # Get variables for negative flow and head difference. + qn, dhn = var(wm, n, :qn_pipe, a), var(wm, n, :dhn_pipe, a) + + # Loop over breakpoints strictly between the lower and upper variable bounds. + for pt in range(0.0, stop = JuMP.upper_bound(qn), length = num_breakpoints+2)[2:end-1] + # Add a linear outer approximation of the convex relaxation at `pt`. + lhs = _get_head_loss_oa_binary(qn, 1.0 - y, pt, exponent) + c = JuMP.@constraint(wm.model, r * lhs <= inv(L) * dhn) + append!(con(wm, n, :pipe_head_loss)[a], [c]) + end + + # Add linear upper bounds on the above outer approximations. + rhs = r * JuMP.upper_bound(qn)^(exponent - 1.0) * qn + c = JuMP.@constraint(wm.model, inv(L) * dhn <= rhs) + + # Append the :pipe_head_loss constraint array. + append!(con(wm, n, :pipe_head_loss)[a], [c]) +end + + +"Pump head gain constraint when the pump status is ambiguous." +function constraint_on_off_pump_head_gain(wm::CQRDWaterModel, n::Int, a::Int, node_fr::Int, node_to::Int, pc::Array{Float64}, q_min_active::Float64) + # Gather pump flow, head gain, and status variables. + qp, g, z = var(wm, n, :qp_pump, a), var(wm, n, :g_pump, a), var(wm, n, :z_pump, a) + + # Define the (relaxed) head gain relationship for the pump. + c_1 = JuMP.@constraint(wm.model, g <= pc[1]*qp^2 + pc[2]*qp + pc[3]*z) + + # Add a linear lower bound on the head gain approximation. + qp_ub = JuMP.upper_bound(qp) + g_1, g_2 = pc[3], pc[1]*qp_ub^2 + pc[2]*qp_ub + pc[3] + g_lb_line = (g_2 - g_1) * inv(qp_ub) * qp + g_1 * z + c_2 = JuMP.@constraint(wm.model, g_lb_line <= g) + + append!(con(wm, n, :on_off_pump_head_gain, a), [c_1, c_2]) +end + + +function objective_wf(wm::CQRDWaterModel) + JuMP.set_objective_sense(wm.model, _MOI.FEASIBILITY_SENSE) +end + + +function objective_owf(wm::CQRDWaterModel) + # Initialize the objective function. + objective = JuMP.AffExpr(0.0) + + for (n, nw_ref) in nws(wm) + # Get common constant parameters. + rho = 1000.0 # Water density (kilogram per cubic meter). + gravity = 9.80665 # Gravitational acceleration (meter per second squared). + constant = rho * gravity * ref(wm, n, :time_step) + + for (a, pump) in nw_ref[:pump] + q_min_active = get(pump, "q_min_active", _q_eps) + + if haskey(pump, "energy_price") + # Get price and pump curve data. + price = pump["energy_price"] + head_curve = ref(wm, n, :pump, a)["head_curve"] + curve_fun = _get_function_from_head_curve(head_curve) + + # Get flow-related variables and data. + z = var(wm, n, :z_pump, a) + qp, g = var(wm, n, :qp_pump, a), var(wm, n, :g_pump, a) + points = collect(range(q_min_active, stop=JuMP.upper_bound(qp), length=50)) + + # Get pump efficiency data. + if haskey(pump, "efficiency_curve") + # Use the efficiency at the midpoint of the pump flow bounds. + eff_curve = pump["efficiency_curve"] + eff = _calc_efficiencies(points, eff_curve) + else + eff = pump["efficiency"] + end + + # Compute discrete costs from existing points. + flows_cubed = _calc_cubic_flow_values(points, curve_fun) + costs = (constant*price) .* inv.(eff) .* flows_cubed + + # Fit a quadratic function to the above discrete costs. + LsqFit.@. func(x, p) = p[1]*x + p[2] + fit = LsqFit.curve_fit(func, points, costs, zeros(length(costs))) + coeffs = LsqFit.coef(fit) + + # Add the cost corresponding to the current pump's operation. + JuMP.add_to_expression!(objective, coeffs[1]*qp + coeffs[2]*z) + else + Memento.error(_LOGGER, "No cost given for pump \"$(pump["name"])\"") + end + end + end + + return JuMP.@objective(wm.model, _MOI.MIN_SENSE, objective) +end diff --git a/src/form/crd.jl b/src/form/crd.jl index 24d4068c..55c22267 100644 --- a/src/form/crd.jl +++ b/src/form/crd.jl @@ -1,132 +1,138 @@ # Define common CRD (continuous or convex relaxation- and direction-based) implementations # of water distribution network constraints, which use directed flow variables. -function constraint_pipe_head_loss_des(wm::CRDWaterModel, n::Int, a::Int, alpha::Float64, node_fr::Int, node_to::Int, L::Float64, pipe_resistances) - # Collect head difference variables. - dhp, dhn = var(wm, n, :dhp_des_pipe, a), var(wm, n, :dhn_des_pipe, a) - for (r_id, r) in enumerate(pipe_resistances) - # Collect directed flow variables. - qp, qn = var(wm, n, :qp_des_pipe, a)[r_id], var(wm, n, :qn_des_pipe, a)[r_id] +########################################## PIPES ########################################## - # Build the relaxed head loss constraints. - c_p = JuMP.@NLconstraint(wm.model, r * head_loss(qp) <= inv(L) * dhp) - c_n = JuMP.@NLconstraint(wm.model, r * head_loss(qn) <= inv(L) * dhn) - # Append the constraint array. - append!(con(wm, n, :head_loss, a), [c_p, c_n]) - end -end +function constraint_pipe_head_loss(wm::CRDWaterModel, n::Int, a::Int, node_fr::Int, node_to::Int, exponent::Float64, L::Float64, r::Float64) + # Gather directed flow and head difference variables. + qp, qn = var(wm, n, :qp_pipe, a), var(wm, n, :qn_pipe, a) + dhp, dhn = var(wm, n, :dhp_pipe, a), var(wm, n, :dhn_pipe, a) -"Pump head gain constraint when the pump status is ambiguous." -function constraint_pump_head_gain(wm::CRDWaterModel, n::Int, a::Int, node_fr::Int, node_to::Int, pc::Array{Float64}) - # Gather flow and head gain variables. - g = var(wm, n, :g, a) - qp, z = var(wm, n, :qp_pump, a), var(wm, n, :z_pump, a) + # Add relaxed constraints for head loss in the positive and negative directions. + c_1 = JuMP.@NLconstraint(wm.model, r * head_loss(qp) <= inv(L) * dhp) + c_2 = JuMP.@NLconstraint(wm.model, r * head_loss(qn) <= inv(L) * dhn) - # Gather head-related variables and data. - h_i, h_j = var(wm, n, :h, node_fr), var(wm, n, :h, node_to) - dhp, dhn = var(wm, n, :dhp_pump, a), var(wm, n, :dhn_pump, a) - dhp_ub, dhn_ub = JuMP.upper_bound(dhp), JuMP.upper_bound(dhn) + # Add linear upper bounds on the above convex relaxations. + rhs_p = r * JuMP.upper_bound(qp)^(exponent - 1.0) * qp + c_3 = JuMP.@constraint(wm.model, inv(L) * dhp <= rhs_p) + rhs_n = r * JuMP.upper_bound(qn)^(exponent - 1.0) * qn + c_4 = JuMP.@constraint(wm.model, inv(L) * dhn <= rhs_n) - # Define the (relaxed) head gain caused by the pump. - g_expr = pc[1]*qp^2 + pc[2]*qp + pc[3]*z - c = JuMP.@constraint(wm.model, g <= g_expr) # Concavified. + # Append the :pipe_head_loss constraint array. + append!(con(wm, n, :pipe_head_loss)[a], [c_1, c_2, c_3, c_4]) +end - # Append the constraint array. - append!(con(wm, n, :head_gain, a), [c]) - # If the number of breakpoints is not positive, no constraints are added. - pump_breakpoints = get(wm.ext, :pump_breakpoints, 0) - if pump_breakpoints <= 0 return end +function constraint_on_off_pipe_head_loss_des(wm::CRDWaterModel, n::Int, a::Int, exponent::Float64, node_fr::Int, node_to::Int, L::Float64, resistances) + # Collect head difference and flow variables. + dhp, dhn = var(wm, n, :dhp_des_pipe, a), var(wm, n, :dhn_des_pipe, a) + qp, qn = var(wm, n, :qp_des_pipe, a), var(wm, n, :qn_des_pipe, a) - # Gather flow, head gain, and convex combination variables. - lambda, x_pw = [var(wm, n, :lambda_pump), var(wm, n, :x_pw_pump)] + for (r_id, r) in enumerate(resistances) + # Build the relaxed head loss constraints. + c_1 = JuMP.@NLconstraint(wm.model, r * head_loss(qp[r_id]) <= inv(L) * dhp) + c_2 = JuMP.@NLconstraint(wm.model, r * head_loss(qn[r_id]) <= inv(L) * dhn) - # Add the required SOS constraints. - c_1 = JuMP.@constraint(wm.model, sum(lambda[a, :]) == z) - c_2 = JuMP.@constraint(wm.model, sum(x_pw[a, :]) == z) - c_3 = JuMP.@constraint(wm.model, lambda[a, 1] <= x_pw[a, 1]) - c_4 = JuMP.@constraint(wm.model, lambda[a, end] <= x_pw[a, end]) + # Append the :on_off_pipe_head_loss_des constraint array. + append!(con(wm, n, :on_off_pipe_head_loss_des, a), [c_1, c_2]) + end - # Add a constraint for the flow piecewise approximation. - qp_ub = JuMP.upper_bound(qp) - breakpoints = range(0.0, stop=qp_ub, length=pump_breakpoints) - qp_lhs = sum(breakpoints[k] * lambda[a, k] for k in 1:pump_breakpoints) - c_5 = JuMP.@constraint(wm.model, qp_lhs == qp) + # Add linear upper bounds for the positive portion of head loss. + qp_ub = JuMP.upper_bound.(qp) + slopes_p = resistances .* qp_ub.^(exponent - 1.0) + c_3 = JuMP.@constraint(wm.model, inv(L)*dhp <= sum(slopes_p .* qp)) - # Append the constraint array. - append!(con(wm, n, :head_gain, a), [c_1, c_2, c_3, c_4, c_5]) + # Add linear upper bounds for the negative portion of head loss. + qn_ub = JuMP.upper_bound.(qn) + slopes_n = resistances .* qn_ub.^(exponent - 1.0) + c_4 = JuMP.@constraint(wm.model, inv(L)*dhn <= sum(slopes_n .* qn)) - for k in 2:pump_breakpoints-1 - # Add the adjacency constraints for piecewise variables. - adjacency = x_pw[a, k-1] + x_pw[a, k] - c_6_k = JuMP.@constraint(wm.model, lambda[a, k] <= adjacency) - append!(con(wm, n, :head_gain, a), [c_6_k]) - end + # Append the :on_off_pipe_head_loss_des constraint array. + append!(con(wm, n, :on_off_pipe_head_loss_des)[a], [c_3, c_4]) end -function constraint_check_valve_head_loss(wm::CRDWaterModel, n::Int, a::Int, node_fr::Int, node_to::Int, L::Float64, r::Float64) - # Gather flow- and check valve-related variables. - qp, qn = var(wm, n, :qp_check_valve, a), var(wm, n, :qn_check_valve, a) - z = var(wm, n, :z_check_valve, a) - # Gather head variables and upper bound data. - h_i, h_j = var(wm, n, :h, node_fr), var(wm, n, :h, node_to) - dhp, dhn = var(wm, n, :dhp_check_valve, a), var(wm, n, :dhn_check_valve, a) - dhp_ub, dhn_ub = JuMP.upper_bound(dhp), JuMP.upper_bound(dhn) +########################################## PUMPS ########################################## - # Add constraints for flow in the positive and negative directions. - lhs = JuMP.@NLexpression(wm.model, r*head_loss(qp) - inv(L)*dhp) - c_p = JuMP.@NLconstraint(wm.model, lhs <= inv(L) * dhp_ub * (1.0 - z)) - # Append the constraint array. - append!(con(wm, n, :head_loss)[a], [c_p]) +"Instantiate variables associated with modeling each pump's head gain." +function variable_pump_head_gain(wm::CRDWaterModel; nw::Int=wm.cnw, bounded::Bool=true, report::Bool=true) + # Initialize variables for total hydraulic head gain from a pump. + g = var(wm, nw)[:g_pump] = JuMP.@variable(wm.model, [a in ids(wm, nw, :pump)], + base_name="$(nw)_g_pump", lower_bound=0.0, # Pump gain is nonnegative. + start=comp_start_value(ref(wm, nw, :pump, a), "g_pump_start")) + + # Initialize an entry to the solution component dictionary for head gains. + report && sol_component_value(wm, nw, :pump, :g, ids(wm, nw, :pump), g) + + # If the number of breakpoints is not positive, return. + pump_breakpoints = get(wm.ext, :pump_breakpoints, 2) + + # Create weights involved in convex combination constraints. + lambda = var(wm, nw)[:lambda_pump] = JuMP.@variable(wm.model, + [a in ids(wm, nw, :pump), k in 1:pump_breakpoints], + base_name="$(nw)_lambda", lower_bound=0.0, upper_bound=1.0, + start=comp_start_value(ref(wm, nw, :pump, a), "lambda_start", k)) + + # Create binary variables involved in convex combination constraints. + x_pw = var(wm, nw)[:x_pw_pump] = JuMP.@variable(wm.model, + [a in ids(wm, nw, :pump), k in 1:pump_breakpoints-1], + base_name="$(nw)_x_pw", binary=true, + start=comp_start_value(ref(wm, nw, :pump, a), "x_pw_start", k)) end -function constraint_shutoff_valve_head_loss(wm::CRDWaterModel, n::Int, a::Int, node_fr::Int, node_to::Int, L::Float64, r::Float64) - # Gather flow- and shutoff valve-related variables. - qp, qn = var(wm, n, :qp_shutoff_valve, a), var(wm, n, :qn_shutoff_valve, a) - y, z = var(wm, n, :y_shutoff_valve, a), var(wm, n, :z_shutoff_valve, a) - # Gather head variables and upper bound data. - dhp, dhn = var(wm, n, :dhp_shutoff_valve, a), var(wm, n, :dhn_shutoff_valve, a) - dhp_ub, dhn_ub = JuMP.upper_bound(dhp), JuMP.upper_bound(dhn) +"Add constraints associated with modeling a pump's head gain." +function constraint_on_off_pump_head_gain(wm::CRDWaterModel, n::Int, a::Int, node_fr::Int, node_to::Int, pc::Array{Float64}, q_min_active::Float64) + # Get the number of breakpoints for the pump. + pump_breakpoints = get(wm.ext, :pump_breakpoints, 2) - # Add constraints for flow in the positive and negative directions. - lhs_p = JuMP.@NLexpression(wm.model, L * r * head_loss(qp) - dhp) - c_1 = JuMP.@NLconstraint(wm.model, lhs_p <= dhp_ub * (1.0 - y)) - c_2 = JuMP.@NLconstraint(wm.model, lhs_p <= dhp_ub * (1.0 - z)) + # Gather pump flow, head gain, and status variables. + qp, g, z = var(wm, n, :qp_pump, a), var(wm, n, :g_pump, a), var(wm, n, :z_pump, a) - lhs_n = JuMP.@NLexpression(wm.model, L * r * head_loss(qn) - dhn) - c_3 = JuMP.@NLconstraint(wm.model, lhs_n <= dhn_ub * y) - c_4 = JuMP.@NLconstraint(wm.model, lhs_n <= dhn_ub * (1.0 - z)) + # Define the (relaxed) head gain relationship for the pump. + c_1 = JuMP.@constraint(wm.model, g <= pc[1]*qp^2 + pc[2]*qp + pc[3]*z) - # Append the constraint array. - append!(con(wm, n, :head_loss)[a], [c_1, c_2, c_3, c_4]) -end + # Gather flow, head gain, and convex combination variables. + lambda, x_pw = var(wm, n, :lambda_pump), var(wm, n, :x_pw_pump) -function constraint_pipe_head_loss(wm::CRDWaterModel, n::Int, a::Int, node_fr::Int, node_to::Int, alpha::Float64, L::Float64, r::Float64) - # Gather common variables. - qp, qn = var(wm, n, :qp_pipe, a), var(wm, n, :qn_pipe, a) - dhp, dhn = var(wm, n, :dhp_pipe, a), var(wm, n, :dhn_pipe, a) + # Add the required SOS constraints. + c_2 = JuMP.@constraint(wm.model, sum(lambda[a, :]) == z) + c_3 = JuMP.@constraint(wm.model, sum(x_pw[a, :]) == z) + c_4 = JuMP.@constraint(wm.model, lambda[a, 1] <= x_pw[a, 1]) + c_5 = JuMP.@constraint(wm.model, lambda[a, end] <= x_pw[a, end]) - # Add constraints for head loss in the positive and negative directions. - c_p = JuMP.@NLconstraint(wm.model, r * head_loss(qp) <= inv(L) * dhp) - c_n = JuMP.@NLconstraint(wm.model, r * head_loss(qn) <= inv(L) * dhn) + # Add a constraint for the flow piecewise approximation. + breakpoints = range(q_min_active, stop=JuMP.upper_bound(qp), length=pump_breakpoints) + qp_lhs = sum(breakpoints[k] * lambda[a, k] for k in 1:pump_breakpoints) + c_6 = JuMP.@constraint(wm.model, qp_lhs == qp) + + # Add a constraint that lower-bounds the head gain variable. + f = (pc[1] .* breakpoints.^2) .+ (pc[2] .* breakpoints) .+ pc[3] + gain_lb_expr = sum(f[k] .* lambda[a, k] for k in 1:pump_breakpoints) + c_7 = JuMP.@constraint(wm.model, gain_lb_expr <= g) # Append the constraint array. - append!(con(wm, n, :head_loss)[a], [c_p, c_n]) -end + append!(con(wm, n, :on_off_pump_head_gain, a), [c_1, c_2, c_3, c_4, c_5, c_6, c_7]) -function objective_wf(wm::CRDWaterModel) - JuMP.set_objective_sense(wm.model, _MOI.FEASIBILITY_SENSE) + for k in 2:pump_breakpoints-1 + # Add the adjacency constraints for piecewise variables. + adjacency = x_pw[a, k-1] + x_pw[a, k] + c_8_k = JuMP.@constraint(wm.model, lambda[a, k] <= adjacency) + append!(con(wm, n, :on_off_pump_head_gain, a), [c_8_k]) + end end + +######################################## OBJECTIVES ######################################## + + +"Instantiate the objective associated with the Optimal Water Flow problem." function objective_owf(wm::CRDWaterModel) - # If the number of breakpoints is not positive, no objective is added. - pump_breakpoints = get(wm.ext, :pump_breakpoints, 0) - if pump_breakpoints <= 0 return end + # Get the number of breakpoints for the pump. + pump_breakpoints = get(wm.ext, :pump_breakpoints, 2) # Initialize the objective function. objective = JuMP.AffExpr(0.0) diff --git a/src/form/directed_flow.jl b/src/form/directed_flow.jl index a5222a2b..10af1ac3 100644 --- a/src/form/directed_flow.jl +++ b/src/form/directed_flow.jl @@ -4,9 +4,8 @@ # qn is nonzero, qp should be zero. -"Initialize variables associated with flow direction. If this variable is equal -to one, the flow direction is from i to j. If it is equal to zero, the flow -direction is from j to i." +"Initialize variables associated with flow direction. If this variable is equal to one, the +flow direction is from i to j. If it is equal to zero, the flow direction is from j to i." function _variable_component_direction( wm::AbstractDirectedModel, component_name::String; nw::Int=wm.cnw, report::Bool=true) # Store the corresponding component symbol. @@ -100,9 +99,9 @@ function _variable_component_flow( end -"Create flow-related variables common to all directed flow models for edge-type components." +"Create flow-related variables common to all directed flow models for node-connecting components." function variable_flow(wm::AbstractDirectedModel; nw::Int=wm.cnw, bounded::Bool=true, report::Bool=true) - for name in ["check_valve", "pipe", "pressure_reducing_valve", "pump", "shutoff_valve"] + for name in ["pipe", "pump", "regulator", "short_pipe", "valve"] # Create directed flow (`qp` and `qn`) variables for each component. _variable_component_flow(wm, name; nw=nw, bounded=bounded, report=report) @@ -161,349 +160,291 @@ function variable_flow_des(wm::AbstractDirectedModel; nw::Int=wm.cnw, bounded::B # Initialize the solution reporting data structures. report && sol_component_value(wm, nw, :des_pipe, :q, ids(wm, nw, :des_pipe), q) - - # Create resistance binary variables. - variable_resistance(wm, nw=nw) end -function variable_pump_operation(wm::AbstractDirectedModel; nw::Int=wm.cnw, report::Bool=true) - # If the number of breakpoints is not positive, return. - pump_breakpoints = get(wm.ext, :pump_breakpoints, 0) +function constraint_pipe_flow(wm::AbstractDirectedModel, n::Int, a::Int) + # Get common flow variables and associated data. + qp, qn, y = var(wm, n, :qp_pipe, a), var(wm, n, :qn_pipe, a), var(wm, n, :y_pipe, a) - if pump_breakpoints > 0 - # Create weights involved in convex combination constraints. - lambda = var(wm, nw)[:lambda_pump] = JuMP.@variable(wm.model, - [a in ids(wm, nw, :pump), k in 1:pump_breakpoints], - base_name="$(nw)_lambda", lower_bound=0.0, upper_bound=1.0, - start=comp_start_value(ref(wm, nw, :pump, a), "lambda_start", k)) + # Constrain directed flow variables based on direction. + qp_ub, qn_ub = JuMP.upper_bound(qp), JuMP.upper_bound(qn) + c_1 = JuMP.@constraint(wm.model, qp <= qp_ub * y) + c_2 = JuMP.@constraint(wm.model, qn <= qn_ub * (1.0 - y)) - # Create binary variables involved in convex combination constraints. - x_pw = var(wm, nw)[:x_pw_pump] = JuMP.@variable(wm.model, - [a in ids(wm, nw, :pump), k in 1:pump_breakpoints-1], - base_name="$(nw)_x_pw", binary=true, - start=comp_start_value(ref(wm, nw, :pump, a), "x_pw_start", k)) - end + # Append the :pipe_flow constraint array. + append!(con(wm, n, :pipe_flow, a), [c_1, c_2]) end -function constraint_check_valve_common(wm::AbstractDirectedModel, n::Int, a::Int, node_fr::Int, node_to::Int) - # Get common flow variables. - qp, qn = var(wm, n, :qp_check_valve, a), var(wm, n, :qn_check_valve, a) - y, z = var(wm, n, :y_check_valve, a), var(wm, n, :z_check_valve, a) - - # If the check valve is closed, flow must be zero. - c_1 = JuMP.@constraint(wm.model, qp <= JuMP.upper_bound(qp) * z) +function constraint_pipe_head(wm::AbstractDirectedModel, n::Int, a::Int, node_fr::Int, node_to::Int) + # Get common flow variables and associated data. + dhp, dhn, y = var(wm, n, :dhp_pipe, a), var(wm, n, :dhn_pipe, a), var(wm, n, :y_pipe, a) - # Get common head variables and associated data. + # Get head variables for from and to nodes. h_i, h_j = var(wm, n, :h, node_fr), var(wm, n, :h, node_to) - dhp, dhn = var(wm, n, :dhp_check_valve, a), var(wm, n, :dhn_check_valve, a) + + # For pipes, the differences must satisfy lower and upper bounds. dhp_ub, dhn_ub = JuMP.upper_bound(dhp), JuMP.upper_bound(dhn) + c_1 = JuMP.@constraint(wm.model, dhp <= dhp_ub * y) + c_2 = JuMP.@constraint(wm.model, dhn <= dhn_ub * (1.0 - y)) + c_3 = JuMP.@constraint(wm.model, dhp - dhn == h_i - h_j) + + # Append the :pipe_head constraint array. + append!(con(wm, n, :pipe_head, a), [c_1, c_2, c_3]) +end - # When the check valve is closed, positive head loss is not possible. - c_2 = JuMP.@constraint(wm.model, dhp <= dhp_ub * z) - # When the check valve is open, negative head loss is not possible. - c_3 = JuMP.@constraint(wm.model, dhn <= dhn_ub * (1.0 - z)) +function constraint_on_off_pipe_flow_des(wm::AbstractDirectedModel, n::Int, a::Int, resistances) + # Get design pipe direction and status variable references. + y, z = var(wm, n, :y_des_pipe, a), var(wm, n, :z_des_pipe, a) - # Constrain head differences based on direction. - c_4 = JuMP.@constraint(wm.model, dhp <= dhp_ub * y) - c_5 = JuMP.@constraint(wm.model, dhn <= dhn_ub * (1.0 - y)) + # Ensure that only one flow can be nonnegative per solution. + c_1 = JuMP.@constraint(wm.model, sum(z) == 1.0) + append!(con(wm, n, :on_off_pipe_flow_des)[a], [c_1]) - # Constrain directed flows based on direction. - c_6 = JuMP.@constraint(wm.model, qp <= JuMP.upper_bound(qp) * y) + for r_id in 1:length(resistances) + # Get directed flow variables and associated data. + qp, qn = var(wm, n, :qp_des_pipe, a)[r_id], var(wm, n, :qn_des_pipe, a)[r_id] + qp_ub, qn_ub = JuMP.upper_bound(qp), JuMP.upper_bound(qn) - # Relate head differences with head variables - c_7 = JuMP.@constraint(wm.model, dhp - dhn == h_i - h_j) + # Constraint the pipes based on direction and construction status. + c_2 = JuMP.@constraint(wm.model, qp <= qp_ub * y) + c_3 = JuMP.@constraint(wm.model, qn <= qn_ub * (1.0 - y)) + c_4 = JuMP.@constraint(wm.model, qp <= qp_ub * z[r_id]) + c_5 = JuMP.@constraint(wm.model, qn <= qn_ub * z[r_id]) - # Append the constraint array. - append!(con(wm, n, :check_valve, a), [c_1, c_2, c_3, c_4, c_5, c_6, c_7]) + # Append the :on_off_pipe_flow_des constraint array. + append!(con(wm, n, :on_off_pipe_flow_des)[a], [c_2, c_3, c_4, c_5]) + end end -function constraint_sv_common(wm::AbstractDirectedModel, n::Int, a::Int, node_fr::Int, node_to::Int) - # Get flow and shutoff valve status variables. - y, z = var(wm, n, :y_shutoff_valve, a), var(wm, n, :z_shutoff_valve, a) - qp, qn = var(wm, n, :qp_shutoff_valve, a), var(wm, n, :qn_shutoff_valve, a) - qp_ub, qn_ub = JuMP.upper_bound(qp), JuMP.upper_bound(qn) +function constraint_on_off_pipe_head_des(wm::AbstractDirectedModel, n::Int, a::Int, node_fr::Int, node_to::Int) + # Get design pipe direction variable reference. + y = var(wm, n, :y_des_pipe, a) - # Get common head variables and associated data. + # Get directed head variables and associated data. h_i, h_j = var(wm, n, :h, node_fr), var(wm, n, :h, node_to) - dhp, dhn = var(wm, n, :dhp_shutoff_valve, a), var(wm, n, :dhn_shutoff_valve, a) + dhp, dhn = var(wm, n, :dhp_des_pipe, a), var(wm, n, :dhn_des_pipe, a) dhp_ub, dhn_ub = JuMP.upper_bound(dhp), JuMP.upper_bound(dhn) - # Constrain head differences based on direction. + # Constrain the pipe heads based on direction variables. c_1 = JuMP.@constraint(wm.model, dhp <= dhp_ub * y) c_2 = JuMP.@constraint(wm.model, dhn <= dhn_ub * (1.0 - y)) - # Constrain directed flows based on direction. - c_3 = JuMP.@constraint(wm.model, qp <= qp_ub * y) - c_4 = JuMP.@constraint(wm.model, qn <= qn_ub * (1.0 - y)) - - # Constrain directed flows based on shutoff valve status. - c_5 = JuMP.@constraint(wm.model, qp <= qp_ub * z) - c_6 = JuMP.@constraint(wm.model, qn <= qn_ub * z) - - # Relate head differences with head variables. - c_7 = JuMP.@constraint(wm.model, dhp - dhn == h_i - h_j) + # Equate head difference variables with heads. + c_3 = JuMP.@constraint(wm.model, dhp - dhn == h_i - h_j) - # Append the constraint array. - append!(con(wm, n, :sv, a), [c_1, c_2, c_3, c_4, c_5, c_6, c_7]) + # Append the :on_off_pipe_flow_des constraint array. + append!(con(wm, n, :on_off_pipe_head_des)[a], [c_1, c_2, c_3]) end -function constraint_pipe_common(wm::AbstractDirectedModel, n::Int, a::Int, node_fr::Int, node_to::Int, alpha::Float64, L::Float64, r::Float64) - # Get common flow variables and associated data. - y = var(wm, n, :y_pipe, a) - qp, qn = var(wm, n, :qp_pipe, a), var(wm, n, :qn_pipe, a) - qp_ub, qn_ub = JuMP.upper_bound(qp), JuMP.upper_bound(qn) +function constraint_on_off_pump_flow(wm::AbstractDirectedModel, n::Int, a::Int, q_min_active::Float64) + # Get pump status variable. + qp, y, z = var(wm, n, :qp_pump, a), var(wm, n, :z_pump, a), var(wm, n, :z_pump, a) - # Constrain directed flow variables based on direction. - c_1 = JuMP.@constraint(wm.model, qp <= qp_ub * y) - c_2 = JuMP.@constraint(wm.model, qn <= qn_ub * (1.0 - y)) - - # Get common head variables and associated data. - dhp, dhn = var(wm, n, :dhp_pipe, a), var(wm, n, :dhn_pipe, a) - h_i, h_j = var(wm, n, :h, node_fr), var(wm, n, :h, node_to) - dhp_ub, dhn_ub = JuMP.upper_bound(dhp), JuMP.upper_bound(dhn) + # If the pump is inactive, flow must be zero. + qp_ub = JuMP.upper_bound(qp) + c_1 = JuMP.@constraint(wm.model, qp >= q_min_active * z) + c_2 = JuMP.@constraint(wm.model, qp <= qp_ub * z) - # Constrain directed head variables based on direction. - c_3 = JuMP.@constraint(wm.model, dhp <= dhp_ub * y) - c_4 = JuMP.@constraint(wm.model, dhn <= dhn_ub * (1.0 - y)) - c_5 = JuMP.@constraint(wm.model, dhp - dhn == h_i - h_j) + # If the pump is on, flow across the pump must be nonnegative. + c_3 = JuMP.@constraint(wm.model, y >= z) # Append the constraint array. - append!(con(wm, n, :pipe)[a], [c_1, c_2, c_3, c_4, c_5]) + append!(con(wm, n, :on_off_pump_flow, a), [c_1, c_2, c_3]) end -function constraint_prv_common(wm::AbstractDirectedModel, n::Int, a::Int, node_fr::Int, node_to::Int, h_prv::Float64) - # Get common flow variables. - qp = var(wm, n, :qp_pressure_reducing_valve, a) - qn = var(wm, n, :qn_pressure_reducing_valve, a) - y = var(wm, n, :y_pressure_reducing_valve, a) - z = var(wm, n, :z_pressure_reducing_valve, a) +function constraint_on_off_pump_head(wm::AbstractDirectedModel, n::Int, a::Int, node_fr::Int, node_to::Int) + # Get head difference variables for the pump. + dhp, dhn = var(wm, n, :dhp_pump, a), var(wm, n, :dhn_pump, a) - # If the pressure reducing valve is open, flow must be appreciably positive. - c_1 = JuMP.@constraint(wm.model, qp <= JuMP.upper_bound(qp) * z) - c_2 = JuMP.@constraint(wm.model, qp <= JuMP.upper_bound(qp) * y) - c_3 = JuMP.@constraint(wm.model, qp >= _q_eps * z) + # Get pump head gain and status variable. + g, z = var(wm, n, :g_pump, a), var(wm, n, :z_pump, a) - # Get common head variables and associated data. - h_i, h_j = var(wm, n, :h, node_fr), var(wm, n, :h, node_to) - dhp = var(wm, n, :dhp_pressure_reducing_valve, a) - dhn = var(wm, n, :dhn_pressure_reducing_valve, a) + # If the pump is off, decouple the head difference relationship. dhp_ub, dhn_ub = JuMP.upper_bound(dhp), JuMP.upper_bound(dhn) + c_1 = JuMP.@constraint(wm.model, dhp <= dhp_ub * (1.0 - z)) + c_2 = JuMP.@constraint(wm.model, dhn <= g + dhn_ub * (1.0 - z)) + c_3 = JuMP.@constraint(wm.model, dhn >= g) - # When the pressure reducing valve is open, the head at node j is predefined. - h_lb, h_ub = JuMP.lower_bound(h_j), JuMP.upper_bound(h_j) - c_4 = JuMP.@constraint(wm.model, h_j >= (1.0 - z) * h_lb + z * h_prv) - c_5 = JuMP.@constraint(wm.model, h_j <= (1.0 - z) * h_ub + z * h_prv) - - # Constrain directed head differences based on status. - c_6 = JuMP.@constraint(wm.model, dhp <= dhp_ub * z) - c_7 = JuMP.@constraint(wm.model, dhn <= dhn_ub * (1.0 - z)) - c_8 = JuMP.@constraint(wm.model, dhp - dhn == h_i - h_j) + # Equate head difference variables with heads. + h_i, h_j = var(wm, n, :h, node_fr), var(wm, n, :h, node_to) + c_4 = JuMP.@constraint(wm.model, dhp - dhn == h_i - h_j) # Append the constraint array. - append!(con(wm, n, :prv, a), [c_1, c_2, c_3, c_4, c_5, c_6, c_7, c_8]) + append!(con(wm, n, :on_off_pump_head, a), [c_1, c_2, c_3, c_4]) end -"Constraint for lower-bounding relaxations of pump head gain." -function constraint_pump_head_gain_lb(wm::AbstractDirectedModel, n::Int, a::Int, node_fr::Int, node_to::Int, pc::Array{Float64}) - # If the number of breakpoints is not positive, no constraints are added. - pump_breakpoints = get(wm.ext, :pump_breakpoints, 0) - if pump_breakpoints <= 0 return end +function constraint_on_off_regulator_flow(wm::AbstractDirectedModel, n::Int, a::Int, q_min_active::Float64) + # Get regulator flow, status, and direction variables. + qp, z = var(wm, n, :qp_regulator, a), var(wm, n, :z_regulator, a) + y = var(wm, n, :y_regulator, a) # Regulator flow direction. - # Gather flow, head gain, and convex combination variables. - qp, g = var(wm, n, :qp_pump, a), var(wm, n, :g, a) - lambda = var(wm, n, :lambda_pump) + # If the regulator is closed, flow must be zero. + qp_lb, qp_ub = max(JuMP.lower_bound(qp), q_min_active), JuMP.upper_bound(qp) + c_1 = JuMP.@constraint(wm.model, qp >= qp_lb * z) + c_2 = JuMP.@constraint(wm.model, qp <= qp_ub * z) + c_3 = JuMP.@constraint(wm.model, qp <= qp_ub * y) - # Add a constraint that lower-bounds the head gain variable. - breakpoints = range(0.0, stop=JuMP.upper_bound(qp), length=pump_breakpoints) - f = (pc[1] .* breakpoints.^2) .+ (pc[2] .* breakpoints) .+ pc[3] - gain_lb = sum(f[k] .* lambda[a, k] for k in 1:pump_breakpoints) - c = JuMP.@constraint(wm.model, gain_lb <= g) - append!(con(wm, n, :head_gain, a), [c]) + # Append the :on_off_regulator_flow constraint array. + append!(con(wm, n, :on_off_regulator_flow, a), [c_1, c_2, c_3]) end -function constraint_pump_common(wm::AbstractDirectedModel, n::Int, a::Int, node_fr::Int, node_to::Int, pc::Array{Float64}) - # Gather variables common to all formulations involving pumps. - y, z = var(wm, n, :y_pump, a), var(wm, n, :z_pump, a) - qp, g = var(wm, n, :qp_pump, a), var(wm, n, :g, a) - dhp, dhn = var(wm, n, :dhp_pump, a), var(wm, n, :dhn_pump, a) - - # If the pump is off, flow across the pump must be zero. - qp_ub = JuMP.upper_bound(qp) - c_1 = JuMP.@constraint(wm.model, qp <= qp_ub * z) - c_2 = JuMP.@constraint(wm.model, qp >= _q_eps * z) +function constraint_on_off_regulator_head(wm::AbstractDirectedModel, n::Int, a::Int, node_fr::Int, node_to::Int, head_setting::Float64) + # Get regulator direction and status variable. + y, z = var(wm, n, :y_regulator, a), var(wm, n, :z_regulator, a) - # If the pump is off, decouple the head difference relationship. - dhp_ub, dhn_ub = JuMP.upper_bound(dhp), JuMP.upper_bound(dhn) - c_3 = JuMP.@constraint(wm.model, dhp <= dhp_ub * (1.0 - z)) - c_4 = JuMP.@constraint(wm.model, dhn <= g + dhn_ub * (1.0 - z)) - c_5 = JuMP.@constraint(wm.model, dhn >= g) + # Get common head variables and associated data. + h_i, h_j = var(wm, n, :h, node_fr), var(wm, n, :h, node_to) + dhp, dhn = var(wm, n, :dhp_regulator, a), var(wm, n, :dhn_regulator, a) - # If the pump is on, flow across the pump must be nonnegative. - c_6 = JuMP.@constraint(wm.model, y >= z) + # When the pressure reducing valve is open, the head at node j is predefined. + h_lb, h_ub = JuMP.lower_bound(h_j), JuMP.upper_bound(h_j) + c_1 = JuMP.@constraint(wm.model, h_j >= (1.0 - z) * h_lb + z * head_setting) + c_2 = JuMP.@constraint(wm.model, h_j <= (1.0 - z) * h_ub + z * head_setting) - # Equate head difference variables with heads. - h_i, h_j = var(wm, n, :h, node_fr), var(wm, n, :h, node_to) + # Constrain directed head differences based on status. + dhp_ub, dhn_ub = JuMP.upper_bound(dhp), JuMP.upper_bound(dhn) + c_3 = JuMP.@constraint(wm.model, dhp <= dhp_ub * z) + c_4 = JuMP.@constraint(wm.model, dhn <= dhn_ub * (1.0 - z)) + c_5 = JuMP.@constraint(wm.model, dhp <= dhp_ub * y) + c_6 = JuMP.@constraint(wm.model, dhn <= dhn_ub * (1.0 - y)) c_7 = JuMP.@constraint(wm.model, dhp - dhn == h_i - h_j) - # Add a linear lower bound on the head gain approximation. - g_1, g_2 = pc[3], pc[1]*qp_ub^2 + pc[2]*qp_ub + pc[3] - g_lb_line = (g_2 - g_1) * inv(qp_ub) * qp + g_1 * z - c_8 = JuMP.@constraint(wm.model, g_lb_line <= g) - - # Append the constraint array. - append!(con(wm, n, :pump, a), [c_1, c_2, c_3, c_4, c_5, c_6, c_7, c_8]) + # Append the :on_off_regulator_head constraint array. + append!(con(wm, n, :on_off_regulator_head, a), [c_1, c_2, c_3, c_4, c_5, c_6, c_7]) end -function constraint_flow_direction_selection_des(wm::AbstractDirectedModel, n::Int, a::Int, pipe_resistances) - y = var(wm, n, :y_des_pipe, a) - - for r_id in 1:length(pipe_resistances) - qp, qn = var(wm, n, :qp_des_pipe, a)[r_id], var(wm, n, :qn_des_pipe, a)[r_id] - c_p = JuMP.@constraint(wm.model, qp <= JuMP.upper_bound(qp) * y) - c_n = JuMP.@constraint(wm.model, qn <= JuMP.upper_bound(qn) * (1.0 - y)) - append!(con(wm, n, :head_loss)[a], [c_p, c_n]) - end -end +function constraint_on_off_valve_flow(wm::AbstractDirectedModel, n::Int, a::Int) + # Get valve flow, direction, and status variables. + qp, qn = var(wm, n, :qp_valve, a), var(wm, n, :qn_valve, a) + y, z = var(wm, n, :y_valve, a), var(wm, n, :z_valve, a) + # The valve flow is constrained by direction and status. + qp_ub, qn_ub = JuMP.upper_bound(qp), JuMP.upper_bound(qn) + c_1 = JuMP.@constraint(wm.model, qp <= qp_ub * y) + c_2 = JuMP.@constraint(wm.model, qn <= qn_ub * (1.0 - y)) + c_3 = JuMP.@constraint(wm.model, qp <= qp_ub * z) + c_4 = JuMP.@constraint(wm.model, qn <= qn_ub * z) -function constraint_head_loss_ub_cv(wm::AbstractDirectedModel, n::Int, a::Int, alpha::Float64, L::Float64, r::Float64) - qp, dhp = var(wm, n, :qp_check_valve, a), var(wm, n, :dhp_check_valve, a) - rhs_p = r * JuMP.upper_bound(qp)^(alpha - 1.0) * qp - c_p = JuMP.@constraint(wm.model, inv(L) * dhp <= rhs_p) - append!(con(wm, n, :head_loss)[a], [c_p]) + # Append the constraint array. + append!(con(wm, n, :on_off_valve_flow, a), [c_1, c_2, c_3, c_4]) end -function constraint_shutoff_valve_head_loss_ub(wm::AbstractDirectedModel, n::Int, a::Int, alpha::Float64, L::Float64, r::Float64) - z = var(wm, n, :z_shutoff_valve, a) - - qp, dhp = var(wm, n, :qp_shutoff_valve, a), var(wm, n, :dhp_shutoff_valve, a) - rhs_p = JuMP.upper_bound(dhp) * (1.0 - z) + L*r*JuMP.upper_bound(qp)^(alpha - 1.0) * qp - c_p = JuMP.@constraint(wm.model, dhp <= rhs_p) - - qn, dhn = var(wm, n, :qn_shutoff_valve, a), var(wm, n, :dhn_shutoff_valve, a) - rhs_n = JuMP.upper_bound(dhn) * (1.0 - z) + L*r*JuMP.upper_bound(qn)^(alpha - 1.0) * qn - c_n = JuMP.@constraint(wm.model, dhn <= rhs_n) - - append!(con(wm, n, :head_loss)[a], [c_p, c_n]) -end - +function constraint_on_off_valve_head(wm::AbstractDirectedModel, n::Int, a::Int, node_fr::Int, node_to::Int) + # Get valve direction and status variables. + y, z = var(wm, n, :y_valve, a), var(wm, n, :z_valve, a) -function constraint_pipe_head_loss_ub(wm::AbstractDirectedModel, n::Int, a::Int, alpha::Float64, L::Float64, r::Float64) - qp, dhp = var(wm, n, :qp_pipe, a), var(wm, n, :dhp_pipe, a) - rhs_p = r * JuMP.upper_bound(qp)^(alpha - 1.0) * qp - c_p = JuMP.@constraint(wm.model, inv(L) * dhp <= rhs_p) + # Get head variables for from and to nodes. + h_i, h_j = var(wm, n, :h, node_fr), var(wm, n, :h, node_to) + dhp, dhn = var(wm, n, :dhp_valve, a), var(wm, n, :dhn_valve, a) - qn, dhn = var(wm, n, :qn_pipe, a), var(wm, n, :dhn_pipe, a) - rhs_n = r * JuMP.upper_bound(qn)^(alpha - 1.0) * qn - c_n = JuMP.@constraint(wm.model, inv(L) * dhn <= rhs_n) + # For valves, head differences must satisfy lower and upper bounds. + dhp_ub, dhn_ub = JuMP.upper_bound(dhp), JuMP.upper_bound(dhn) + c_1 = JuMP.@constraint(wm.model, dhp <= dhp_ub * y) + c_2 = JuMP.@constraint(wm.model, dhn <= dhn_ub * (1.0 - y)) + c_3 = JuMP.@constraint(wm.model, dhp <= dhp_ub * (1.0 - z)) + c_4 = JuMP.@constraint(wm.model, dhn <= dhn_ub * (1.0 - z)) + c_5 = JuMP.@constraint(wm.model, dhp - dhn == h_i - h_j) - append!(con(wm, n, :head_loss)[a], [c_p, c_n]) + # Append the constraint array. + append!(con(wm, n, :on_off_valve_head, a), [c_1, c_2, c_3, c_4, c_5]) end -function constraint_pipe_head_loss_ub_des(wm::AbstractDirectedModel, n::Int, a::Int, alpha::Float64, L::Float64, pipe_resistances) - dhp = var(wm, n, :dhp_des_pipe, a) - qp_des_pipe = var(wm, n, :qp_des_pipe, a) - qp_des_pipe_ub = JuMP.upper_bound.(qp_des_pipe) - slopes_p = pipe_resistances .* qp_des_pipe_ub.^(alpha - 1.0) - c_p = JuMP.@constraint(wm.model, inv(L)*dhp <= sum(slopes_p .* qp_des_pipe)) +function constraint_short_pipe_flow(wm::AbstractDirectedModel, n::Int, a::Int) + # Get short pipe flow and direction variables. + qp, qn = var(wm, n, :qp_short_pipe, a), var(wm, n, :qn_short_pipe, a) + y = var(wm, n, :y_short_pipe, a) # Binary direction variable. - dhn = var(wm, n, :dhn_des_pipe, a) - qn_des_pipe = var(wm, n, :qn_des_pipe, a) - qn_des_pipe_ub = JuMP.upper_bound.(qn_des_pipe) - slopes_n = pipe_resistances .* qn_des_pipe_ub.^(alpha - 1.0) - c_n = JuMP.@constraint(wm.model, inv(L)*dhn <= sum(slopes_n .* qn_des_pipe)) + # The valve flow is constrained by direction and status. + qp_ub, qn_ub = JuMP.upper_bound(qp), JuMP.upper_bound(qn) + c_1 = JuMP.@constraint(wm.model, qp <= qp_ub * y) + c_2 = JuMP.@constraint(wm.model, qn <= qn_ub * (1.0 - y)) - append!(con(wm, n, :head_loss)[a], [c_p, c_n]) + # Append the :short_pipe_flow constraint array. + append!(con(wm, n, :short_pipe_flow, a), [c_1, c_2]) end -function constraint_resistance_selection_des(wm::AbstractDirectedModel, n::Int, a::Int, pipe_resistances) - c = JuMP.@constraint(wm.model, sum(var(wm, n, :x_res, a)) == 1.0) - append!(con(wm, n, :head_loss)[a], [c]) - - for (r_id, r) in enumerate(pipe_resistances) - x_res = var(wm, n, :x_res, a)[r_id] +function constraint_short_pipe_head(wm::AbstractDirectedModel, n::Int, a::Int, node_fr::Int, node_to::Int) + # Get short pipe direction variable. + y = var(wm, n, :y_short_pipe, a) - qp_des_pipe = var(wm, n, :qp_des_pipe, a)[r_id] - qp_ub = JuMP.upper_bound(qp_des_pipe) - c_p = JuMP.@constraint(wm.model, qp_des_pipe <= qp_ub * x_res) + # Get head variables for from and to nodes. + h_i, h_j = var(wm, n, :h, node_fr), var(wm, n, :h, node_to) + dhp, dhn = var(wm, n, :dhp_short_pipe, a), var(wm, n, :dhn_short_pipe, a) - qn_des_pipe = var(wm, n, :qn_des_pipe, a)[r_id] - qn_ub = JuMP.upper_bound(qn_des_pipe) - c_n = JuMP.@constraint(wm.model, qn_des_pipe <= qn_ub * x_res) + # For short pipes, head differences must satisfy lower and upper bounds. + dhp_ub, dhn_ub = JuMP.upper_bound(dhp), JuMP.upper_bound(dhn) + c_1 = JuMP.@constraint(wm.model, h_i - h_j == 0.0) + c_2 = JuMP.@constraint(wm.model, dhp <= dhp_ub * y) + c_3 = JuMP.@constraint(wm.model, dhn <= dhn_ub * (1.0 - y)) + c_4 = JuMP.@constraint(wm.model, dhp - dhn == h_i - h_j) - append!(con(wm, n, :head_loss)[a], [c_p, c_n]) - end + # Append the :short_pipe_head constraint array. + append!(con(wm, n, :short_pipe_head, a), [c_1, c_2, c_3, c_4]) end function _gather_directionality_data( - wm::AbstractDirectedModel, n::Int, check_valve_fr::Array{Int64,1}, - check_valve_to::Array{Int64,1}, pipe_fr::Array{Int64,1}, pipe_to::Array{Int64,1}, + wm::AbstractDirectedModel, n::Int, pipe_fr::Array{Int64,1}, pipe_to::Array{Int64,1}, des_pipe_fr::Array{Int64,1}, des_pipe_to::Array{Int64,1}, pump_fr::Array{Int64,1}, - pump_to::Array{Int64,1}, pressure_reducing_valve_fr::Array{Int64,1}, - pressure_reducing_valve_to::Array{Int64,1}, shutoff_valve_fr::Array{Int64,1}, - shutoff_valve_to::Array{Int64,1}) + pump_to::Array{Int64,1}, regulator_fr::Array{Int64,1}, regulator_to::Array{Int64,1}, + short_pipe_fr::Array{Int64,1}, short_pipe_to::Array{Int64,1}, valve_fr::Array{Int64,1}, + valve_to::Array{Int64,1}) # Collect direction variable references per component. - y_check_valve = var(wm, n, :y_check_valve) - y_pipe = var(wm, n, :y_pipe) - y_des_pipe = var(wm, n, :y_des_pipe) - y_pump = var(wm, n, :y_pump) - y_pressure_reducing_valve = var(wm, n, :y_pressure_reducing_valve) - y_shutoff_valve = var(wm, n, :y_shutoff_valve) - - sum_out = JuMP.@expression(wm.model, - sum(y_check_valve[a] for a in check_valve_fr) + - sum(y_pipe[a] for a in pipe_fr) + - sum(y_des_pipe[a] for a in des_pipe_fr) + - sum(y_pump[a] for a in pump_fr) + - sum(y_pressure_reducing_valve[a] for a in pressure_reducing_valve_fr) + - sum(y_shutoff_valve[a] for a in shutoff_valve_fr)) + y_pipe, y_des_pipe = var(wm, n, :y_pipe), var(wm, n, :y_des_pipe) + y_pump, y_regulator = var(wm, n, :y_pump), var(wm, n, :y_regulator) + y_short_pipe, y_valve = var(wm, n, :y_short_pipe), var(wm, n, :y_valve) sum_in = JuMP.@expression(wm.model, - sum(y_check_valve[a] for a in check_valve_to) + sum(y_pipe[a] for a in pipe_to) + sum(y_des_pipe[a] for a in des_pipe_to) + sum(y_pump[a] for a in pump_to) + - sum(y_pressure_reducing_valve[a] for a in pressure_reducing_valve_to) + - sum(y_shutoff_valve[a] for a in shutoff_valve_to)) + sum(y_regulator[a] for a in regulator_to) + + sum(y_short_pipe[a] for a in short_pipe_to) + + sum(y_valve[a] for a in valve_to)) + + sum_out = JuMP.@expression(wm.model, + sum(y_pipe[a] for a in pipe_fr) + + sum(y_des_pipe[a] for a in des_pipe_fr) + + sum(y_pump[a] for a in pump_fr) + + sum(y_regulator[a] for a in regulator_fr) + + sum(y_short_pipe[a] for a in short_pipe_fr) + + sum(y_valve[a] for a in valve_fr)) # Get the in degree of node `i`. - in_length = length(check_valve_to) + length(pipe_to) + length(des_pipe_to) + - length(pump_to) + length(pressure_reducing_valve_to) + length(shutoff_valve_to) + in_length = length(pipe_to) + length(des_pipe_to) + length(pump_to) + + length(regulator_to) + length(short_pipe_to) + length(valve_to) # Get the out degree of node `i`. - out_length = length(check_valve_fr) + length(pipe_fr) + length(des_pipe_fr) + - length(pump_fr) + length(pressure_reducing_valve_fr) + length(shutoff_valve_fr) + out_length = length(pipe_fr) + length(des_pipe_fr) + length(pump_fr) + + length(regulator_fr) + length(short_pipe_fr) + length(valve_fr) return sum_in, sum_out, in_length, out_length end "Constraint to ensure at least one direction is set to take flow away from a source." -function constraint_node_directionality( - wm::AbstractDirectedModel, n::Int, i::Int, check_valve_fr::Array{Int64,1}, - check_valve_to::Array{Int64,1}, pipe_fr::Array{Int64,1}, pipe_to::Array{Int64,1}, - des_pipe_fr::Array{Int64,1}, des_pipe_to::Array{Int64,1}, pump_fr::Array{Int64,1}, - pump_to::Array{Int64,1}, pressure_reducing_valve_fr::Array{Int64,1}, - pressure_reducing_valve_to::Array{Int64,1}, shutoff_valve_fr::Array{Int64,1}, - shutoff_valve_to::Array{Int64,1}) +function constraint_intermediate_directionality( + wm::AbstractDirectedModel, n::Int, i::Int, pipe_fr::Array{Int64,1}, + pipe_to::Array{Int64,1}, des_pipe_fr::Array{Int64,1}, des_pipe_to::Array{Int64,1}, + pump_fr::Array{Int64,1}, pump_to::Array{Int64,1}, regulator_fr::Array{Int64,1}, + regulator_to::Array{Int64,1}, short_pipe_fr::Array{Int64,1}, + short_pipe_to::Array{Int64,1}, valve_fr::Array{Int64,1}, valve_to::Array{Int64,1}) # Gather data required to build the constraint. sum_in, sum_out, in_length, out_length = _gather_directionality_data( - wm, n, check_valve_fr, check_valve_to, pipe_fr, pipe_to, des_pipe_fr, des_pipe_to, - pump_fr, pump_to, pressure_reducing_valve_fr, pressure_reducing_valve_to, - shutoff_valve_fr, shutoff_valve_to) + wm, n, pipe_fr, pipe_to, des_pipe_fr, des_pipe_to, pump_fr, pump_to, regulator_fr, + regulator_to, short_pipe_fr, short_pipe_to, valve_fr, valve_to) # Add the node directionality constraint. if out_length == 1 && in_length == 1 @@ -518,39 +459,35 @@ end "Constraint to ensure at least one direction is set to take flow away from a source." function constraint_source_directionality( - wm::AbstractDirectedModel, n::Int, i::Int, check_valve_fr::Array{Int64,1}, - check_valve_to::Array{Int64,1}, pipe_fr::Array{Int64,1}, pipe_to::Array{Int64,1}, - des_pipe_fr::Array{Int64,1}, des_pipe_to::Array{Int64,1}, pump_fr::Array{Int64,1}, - pump_to::Array{Int64,1}, pressure_reducing_valve_fr::Array{Int64,1}, - pressure_reducing_valve_to::Array{Int64,1}, shutoff_valve_fr::Array{Int64,1}, - shutoff_valve_to::Array{Int64,1}) + wm::AbstractDirectedModel, n::Int, i::Int, pipe_fr::Array{Int64,1}, + pipe_to::Array{Int64,1}, des_pipe_fr::Array{Int64,1}, des_pipe_to::Array{Int64,1}, + pump_fr::Array{Int64,1}, pump_to::Array{Int64,1}, regulator_fr::Array{Int64,1}, + regulator_to::Array{Int64,1}, short_pipe_fr::Array{Int64,1}, + short_pipe_to::Array{Int64,1}, valve_fr::Array{Int64,1}, valve_to::Array{Int64,1}) # Gather data required to build the constraint. sum_in, sum_out, in_length, out_length = _gather_directionality_data( - wm, n, check_valve_fr, check_valve_to, pipe_fr, pipe_to, des_pipe_fr, des_pipe_to, - pump_fr, pump_to, pressure_reducing_valve_fr, pressure_reducing_valve_to, - shutoff_valve_fr, shutoff_valve_to) + wm, n, pipe_fr, pipe_to, des_pipe_fr, des_pipe_to, pump_fr, pump_to, regulator_fr, + regulator_to, short_pipe_fr, short_pipe_to, valve_fr, valve_to) # Add the source flow direction constraint. c = JuMP.@constraint(wm.model, sum_out - sum_in >= 1.0 - in_length) - con(wm, n, :source_directionality)[i] = c + con(wm, n, :node_directionality)[i] = c end "Constraint to ensure at least one direction is set to take flow to a junction with demand." function constraint_sink_directionality( - wm::AbstractDirectedModel, n::Int, i::Int, check_valve_fr::Array{Int64,1}, - check_valve_to::Array{Int64,1}, pipe_fr::Array{Int64,1}, pipe_to::Array{Int64,1}, - des_pipe_fr::Array{Int64,1}, des_pipe_to::Array{Int64,1}, pump_fr::Array{Int64,1}, - pump_to::Array{Int64,1}, pressure_reducing_valve_fr::Array{Int64,1}, - pressure_reducing_valve_to::Array{Int64,1}, shutoff_valve_fr::Array{Int64,1}, - shutoff_valve_to::Array{Int64,1}) + wm::AbstractDirectedModel, n::Int, i::Int, pipe_fr::Array{Int64,1}, + pipe_to::Array{Int64,1}, des_pipe_fr::Array{Int64,1}, des_pipe_to::Array{Int64,1}, + pump_fr::Array{Int64,1}, pump_to::Array{Int64,1}, regulator_fr::Array{Int64,1}, + regulator_to::Array{Int64,1}, short_pipe_fr::Array{Int64,1}, + short_pipe_to::Array{Int64,1}, valve_fr::Array{Int64,1}, valve_to::Array{Int64,1}) # Gather data required to build the constraint. sum_in, sum_out, in_length, out_length = _gather_directionality_data( - wm, n, check_valve_fr, check_valve_to, pipe_fr, pipe_to, des_pipe_fr, des_pipe_to, - pump_fr, pump_to, pressure_reducing_valve_fr, pressure_reducing_valve_to, - shutoff_valve_fr, shutoff_valve_to) + wm, n, pipe_fr, pipe_to, des_pipe_fr, des_pipe_to, pump_fr, pump_to, regulator_fr, + regulator_to, short_pipe_fr, short_pipe_to, valve_fr, valve_to) # Add the sink flow direction constraint. c = JuMP.@constraint(wm.model, sum_in - sum_out >= 1.0 - out_length) - con(wm, n, :sink_directionality)[i] = c + con(wm, n, :node_directionality)[i] = c end diff --git a/src/form/la.jl b/src/form/la.jl index 3b687515..be8f0a7f 100644 --- a/src/form/la.jl +++ b/src/form/la.jl @@ -1,87 +1,60 @@ # Define LA (linear approximation-based) implementations of water models. function variable_flow_piecewise_weights(wm::LAWaterModel; nw::Int=wm.cnw, report::Bool=false) - pipe_breakpoints = get(wm.ext, :pipe_breakpoints, 0) - - if pipe_breakpoints > 0 - # Create weights involved in convex combination constraints for pipes. - lambda_pipe = var(wm, nw)[:lambda_pipe] = JuMP.@variable(wm.model, - [a in ids(wm, nw, :pipe), k in 1:pipe_breakpoints], - base_name="$(nw)_lambda", lower_bound=0.0, upper_bound=1.0, - start=comp_start_value(ref(wm, nw, :pipe, a), "lambda_start", k)) - - # Create weights involved in convex combination constraints for check valves. - lambda_check_valve = var(wm, nw)[:lambda_check_valve] = JuMP.@variable(wm.model, - [a in ids(wm, nw, :check_valve), k in 1:pipe_breakpoints], - base_name="$(nw)_lambda", lower_bound=0.0, upper_bound=1.0, - start=comp_start_value(ref(wm, nw, :check_valve, a), "lambda_start", k)) - - # Create weights involved in convex combination constraints for shutoff valves. - lambda_shutoff_valve = var(wm, nw)[:lambda_shutoff_valve] = JuMP.@variable(wm.model, - [a in ids(wm, nw, :shutoff_valve), k in 1:pipe_breakpoints], - base_name="$(nw)_lambda", lower_bound=0.0, upper_bound=1.0, - start=comp_start_value(ref(wm, nw, :shutoff_valve, a), "lambda_start", k)) - - # Create weights involved in convex combination constraints. - n_r = Dict(a=>length(ref(wm, nw, :resistance, a)) for a in ids(wm, nw, :des_pipe)) - - lambda_des_pipe = var(wm, nw)[:lambda_des_pipe] = JuMP.@variable(wm.model, - [a in ids(wm, nw, :des_pipe), r in 1:n_r[a], k in 1:pipe_breakpoints], - base_name="$(nw)_lambda", lower_bound=0.0, upper_bound=1.0, - start=comp_start_value(ref(wm, nw, :des_pipe, a), "lambda_start", r)) - end - - pump_breakpoints = get(wm.ext, :pump_breakpoints, 0) - - if pump_breakpoints > 0 - # Create weights involved in convex combination constraints for pumps. - lambda_pump = var(wm, nw)[:lambda_pump] = JuMP.@variable(wm.model, - [a in ids(wm, nw, :pump), k in 1:pump_breakpoints], - base_name="$(nw)_lambda", lower_bound=0.0, upper_bound=1.0, - start=comp_start_value(ref(wm, nw, :pump, a), "lambda_start", k)) - end + # Get the number of breakpoints for the pipe. + pipe_breakpoints = get(wm.ext, :pipe_breakpoints, 2) + + # Create weights involved in convex combination constraints for pipes. + lambda_pipe = var(wm, nw)[:lambda_pipe] = JuMP.@variable(wm.model, + [a in ids(wm, nw, :pipe), k in 1:pipe_breakpoints], + base_name="$(nw)_lambda", lower_bound=0.0, upper_bound=1.0, + start=comp_start_value(ref(wm, nw, :pipe, a), "lambda_start", k)) + + # Create weights involved in convex combination constraints. + n_r = Dict(a=>length(ref(wm, nw, :resistance, a)) for a in ids(wm, nw, :des_pipe)) + lambda_des_pipe = var(wm, nw)[:lambda_des_pipe] = JuMP.@variable(wm.model, + [a in ids(wm, nw, :des_pipe), r in 1:n_r[a], k in 1:pipe_breakpoints], + base_name="$(nw)_lambda", lower_bound=0.0, upper_bound=1.0, + start=comp_start_value(ref(wm, nw, :des_pipe, a), "lambda_start", r)) + + # Get the number of breakpoints for the pump. + pump_breakpoints = get(wm.ext, :pump_breakpoints, 2) + + # Create weights involved in convex combination constraints for pumps. + lambda_pump = var(wm, nw)[:lambda_pump] = JuMP.@variable(wm.model, + [a in ids(wm, nw, :pump), k in 1:pump_breakpoints], + base_name="$(nw)_lambda", lower_bound=0.0, upper_bound=1.0, + start=comp_start_value(ref(wm, nw, :pump, a), "lambda_start", k)) end function variable_flow_piecewise_adjacency(wm::LAWaterModel; nw::Int=wm.cnw, report::Bool=false) - pipe_breakpoints = get(wm.ext, :pipe_breakpoints, 0) - - if pipe_breakpoints > 0 - # Create binary variables for pipe convex combination constraints. - x_pw_pipe = var(wm, nw)[:x_pw_pipe] = JuMP.@variable(wm.model, - [a in ids(wm, nw, :pipe), k in 1:pipe_breakpoints-1], base_name="$(nw)_x_pw", binary=true, - start=comp_start_value(ref(wm, nw, :pipe, a), "x_pw_start")) - - # Create binary variables for check valve convex combination constraints. - x_pw_check_valve = var(wm, nw)[:x_pw_check_valve] = JuMP.@variable(wm.model, - [a in ids(wm, nw, :check_valve), k in 1:pipe_breakpoints-1], base_name="$(nw)_x_pw", binary=true, - start=comp_start_value(ref(wm, nw, :check_valve, a), "x_pw_start")) - - # Create binary variables for shutoff valve convex combination constraints. - x_pw_shutoff_valve = var(wm, nw)[:x_pw_shutoff_valve] = JuMP.@variable(wm.model, - [a in ids(wm, nw, :shutoff_valve), k in 1:pipe_breakpoints-1], base_name="$(nw)_x_pw", binary=true, - start=comp_start_value(ref(wm, nw, :shutoff_valve, a), "x_pw_start")) - - # Create binary variables for design pipe convex combination constraints. - x_pw_des_pipe = var(wm, nw)[:x_pw_des_pipe] = JuMP.@variable(wm.model, - [a in ids(wm, nw, :des_pipe), k in 1:pipe_breakpoints-1], base_name="$(nw)_x_pw", binary=true, - start=comp_start_value(ref(wm, nw, :des_pipe, a), "x_pw_start")) - end - - pump_breakpoints = get(wm.ext, :pump_breakpoints, 0) - - if pump_breakpoints > 0 - # Create binary variables involved in convex combination constraints for pumps. - x_pw_pump = var(wm, nw)[:x_pw_pump] = JuMP.@variable(wm.model, - [a in ids(wm, nw, :pump), k in 1:pump_breakpoints-1], base_name="$(nw)_x_pw", - binary=true, start=comp_start_value(ref(wm, nw, :pump, a), "x_pw_start")) - end + # Get the number of breakpoints for the pipe. + pipe_breakpoints = get(wm.ext, :pipe_breakpoints, 2) + + # Create binary variables for pipe convex combination constraints. + x_pw_pipe = var(wm, nw)[:x_pw_pipe] = JuMP.@variable(wm.model, + [a in ids(wm, nw, :pipe), k in 1:pipe_breakpoints-1], base_name="$(nw)_x_pw", binary=true, + start=comp_start_value(ref(wm, nw, :pipe, a), "x_pw_start")) + + # Create binary variables for design pipe convex combination constraints. + x_pw_des_pipe = var(wm, nw)[:x_pw_des_pipe] = JuMP.@variable(wm.model, + [a in ids(wm, nw, :des_pipe), k in 1:pipe_breakpoints-1], base_name="$(nw)_x_pw", binary=true, + start=comp_start_value(ref(wm, nw, :des_pipe, a), "x_pw_start")) + + # Get the number of breakpoints for the pump. + pump_breakpoints = get(wm.ext, :pump_breakpoints, 2) + + # Create binary variables involved in convex combination constraints for pumps. + x_pw_pump = var(wm, nw)[:x_pw_pump] = JuMP.@variable(wm.model, + [a in ids(wm, nw, :pump), k in 1:pump_breakpoints-1], base_name="$(nw)_x_pw", + binary=true, start=comp_start_value(ref(wm, nw, :pump, a), "x_pw_start")) end "Creates flow variables for `LA` formulations (`q`, `lambda`, `x_pw`)." function variable_flow(wm::LAWaterModel; nw::Int=wm.cnw, bounded::Bool=true, report::Bool=true) - for name in ["check_valve", "pipe", "pressure_reducing_valve", "pump", "shutoff_valve"] + for name in ["pipe", "pump", "regulator", "short_pipe", "valve"] # Create directed flow (`qp` and `qn`) variables for each component. variable_component_flow(wm, name; nw=nw, bounded=bounded, report=report) end @@ -95,16 +68,14 @@ function variable_flow(wm::LAWaterModel; nw::Int=wm.cnw, bounded::Bool=true, rep end -"Adds head loss constraints for check valves in `LA` formulations." -function constraint_check_valve_head_loss(wm::LAWaterModel, n::Int, a::Int, node_fr::Int, node_to::Int, L::Float64, r::Float64) - # If the number of breakpoints is not positive, no constraints are added. - pipe_breakpoints = get(wm.ext, :pipe_breakpoints, 0) - if pipe_breakpoints <= 0 return end +"Pump head gain constraint when the pump status is ambiguous." +function constraint_on_off_pump_head_gain(wm::LAWaterModel, n::Int, a::Int, node_fr::Int, node_to::Int, curve_fun::Array{Float64}, q_min_active::Float64) + # Get the number of breakpoints for the pump. + pump_breakpoints = get(wm.ext, :pump_breakpoints, 2) # Gather common variables. - q, z = var(wm, n, :q_check_valve, a), var(wm, n, :z_check_valve, a) - lambda, x_pw = var(wm, n, :lambda_check_valve), var(wm, n, :x_pw_check_valve) - h_i, h_j = var(wm, n, :h, node_fr), var(wm, n, :h, node_to) + q, g, z = var(wm, n, :q_pump, a), var(wm, n, :g_pump, a), var(wm, n, :z_pump, a) + lambda, x_pw = var(wm, n, :lambda_pump), var(wm, n, :x_pw_pump) # Add the required SOS constraints. c_1 = JuMP.@constraint(wm.model, sum(lambda[a, :]) == z) @@ -112,52 +83,45 @@ function constraint_check_valve_head_loss(wm::LAWaterModel, n::Int, a::Int, node c_3 = JuMP.@constraint(wm.model, lambda[a, 1] <= x_pw[a, 1]) c_4 = JuMP.@constraint(wm.model, lambda[a, end] <= x_pw[a, end]) - # Generate a set of uniform flow and head loss breakpoints. - q_lb, q_ub = JuMP.lower_bound(q), JuMP.upper_bound(q) - breakpoints = range(q_lb, stop=q_ub, length=pipe_breakpoints) - f = _calc_head_loss_values(collect(breakpoints), ref(wm, :alpha)) + # Generate a set of uniform flow breakpoints. + breakpoints = range(q_min_active, stop=JuMP.upper_bound(q), length=pump_breakpoints) # Add a constraint for the flow piecewise approximation. - q_lhs = sum(breakpoints[k] * lambda[a, k] for k in 1:pipe_breakpoints) + q_lhs = sum(breakpoints[k] * lambda[a, k] for k in 1:pump_breakpoints) c_5 = JuMP.@constraint(wm.model, q_lhs == q) # Add a constraint for the head gain piecewise approximation. - loss = r .* sum(f[k] .* lambda[a, k] for k in 1:pipe_breakpoints) - lhs = inv(L) * (h_i - h_j) - loss - c_6 = JuMP.@constraint(wm.model, lhs <= 0.0) - dh_lb = JuMP.lower_bound(h_i) - JuMP.upper_bound(h_j) - c_7 = JuMP.@constraint(wm.model, lhs >= inv(L) * (1.0 - z) * dh_lb) + f = _calc_pump_gain_values(collect(breakpoints), curve_fun) + g_lhs = sum(f[k] * lambda[a, k] for k in 1:pump_breakpoints) + c_6 = JuMP.@constraint(wm.model, g_lhs == g) # Append the constraint array with the above-generated constraints. - append!(con(wm, n, :head_loss, a), [c_1, c_2, c_3, c_4, c_5, c_6, c_7]) + append!(con(wm, n, :on_off_pump_head_gain, a), [c_1, c_2, c_3, c_4, c_5, c_6]) - # Add the adjacency constraints. - for k in 2:pipe_breakpoints-1 + for k in 2:pump_breakpoints-1 + # Add adjacency constraints for each interval. adjacency = x_pw[a, k-1] + x_pw[a, k] - c_8_k = JuMP.@constraint(wm.model, lambda[a, k] <= adjacency) - append!(con(wm, n, :head_loss, a), [c_8_k]) + c_7_k = JuMP.@constraint(wm.model, lambda[a, k] <= adjacency) + append!(con(wm, n, :on_off_pump_head_gain, a), [c_7_k]) end end -"Adds head loss constraints for shutoff valves in `LA` formulations." -function constraint_shutoff_valve_head_loss(wm::LAWaterModel, n::Int, a::Int, node_fr::Int, node_to::Int, L::Float64, r::Float64) +function constraint_pipe_head_loss(wm::LAWaterModel, n::Int, a::Int, node_fr::Int, node_to::Int, exponent::Float64, L::Float64, r::Float64) # If the number of breakpoints is not positive, no constraints are added. - pipe_breakpoints = get(wm.ext, :pipe_breakpoints, 0) - if pipe_breakpoints <= 0 return end + pipe_breakpoints = get(wm.ext, :pipe_breakpoints, 2) - # Gather common variables. - q, z = var(wm, n, :q_shutoff_valve, a), var(wm, n, :z_shutoff_valve, a) - lambda, x_pw = var(wm, n, :lambda_shutoff_valve), var(wm, n, :x_pw_shutoff_valve) - h_i, h_j = var(wm, n, :h, node_fr), var(wm, n, :h, node_to) + # Get required variables. + q, h_i, h_j = var(wm, n, :q_pipe, a), var(wm, n, :h, node_fr), var(wm, n, :h, node_to) + lambda, x_pw = var(wm, n, :lambda_pipe), var(wm, n, :x_pw_pipe) # Add the required SOS constraints. - c_1 = JuMP.@constraint(wm.model, sum(lambda[a, :]) == z) - c_2 = JuMP.@constraint(wm.model, sum(x_pw[a, :]) == z) + c_1 = JuMP.@constraint(wm.model, sum(lambda[a, :]) == 1.0) + c_2 = JuMP.@constraint(wm.model, sum(x_pw[a, :]) == 1.0) c_3 = JuMP.@constraint(wm.model, lambda[a, 1] <= x_pw[a, 1]) c_4 = JuMP.@constraint(wm.model, lambda[a, end] <= x_pw[a, end]) - # Generate a set of uniform flow and head loss breakpoints. + # Generate a set of uniform flow breakpoints. q_lb, q_ub = JuMP.lower_bound(q), JuMP.upper_bound(q) breakpoints = range(q_lb, stop=q_ub, length=pipe_breakpoints) @@ -165,87 +129,42 @@ function constraint_shutoff_valve_head_loss(wm::LAWaterModel, n::Int, a::Int, no q_lhs = sum(breakpoints[k] * lambda[a, k] for k in 1:pipe_breakpoints) c_5 = JuMP.@constraint(wm.model, q_lhs == q) - # Add a constraint for the head gain piecewise approximation. - f = _calc_head_loss_values(collect(breakpoints), ref(wm, :alpha)) - lhs = (h_i - h_j) - L * r * sum(f[k] .* lambda[a, k] for k in 1:pipe_breakpoints) - - dh_ub = JuMP.upper_bound(h_i) - JuMP.lower_bound(h_j) - c_6 = JuMP.@constraint(wm.model, lhs <= (1.0 - z) * dh_ub) - dh_lb = JuMP.lower_bound(h_i) - JuMP.upper_bound(h_j) - c_7 = JuMP.@constraint(wm.model, lhs >= (1.0 - z) * dh_lb) + # Add a constraint for the head loss piecewise approximation. + f = _calc_head_loss_values(collect(breakpoints), exponent) + lhs = r * sum(f[k] * lambda[a, k] for k in 1:pipe_breakpoints) + c_6 = JuMP.@constraint(wm.model, lhs == inv(L) * (h_i - h_j)) # Append the constraint array with the above-generated constraints. - append!(con(wm, n, :head_loss, a), [c_1, c_2, c_3, c_4, c_5, c_6, c_7]) + append!(con(wm, n, :pipe_head_loss, a), [c_1, c_2, c_3, c_4, c_5, c_6]) # Add the adjacency constraints. for k in 2:pipe_breakpoints-1 - adjacency = x_pw[a, k-1] + x_pw[a, k] - c_8_k = JuMP.@constraint(wm.model, lambda[a, k] <= adjacency) - append!(con(wm, n, :head_loss, a), [c_8_k]) - end -end - - -"Pump head gain constraint when the pump status is ambiguous." -function constraint_pump_head_gain(wm::LAWaterModel, n::Int, a::Int, node_fr::Int, node_to::Int, curve_fun::Array{Float64}) - # If the number of breakpoints is not positive, no constraints are added. - pump_breakpoints = get(wm.ext, :pump_breakpoints, 0) - if pump_breakpoints <= 0 return end - - # Gather common variables. - q, g, z = var(wm, n, :q_pump, a), var(wm, n, :g, a), var(wm, n, :z_pump, a) - lambda, x_pw = var(wm, n, :lambda_pump), var(wm, n, :x_pw_pump) - - # Add the required SOS constraints. - c_1 = JuMP.@constraint(wm.model, sum(lambda[a, :]) == z) - c_2 = JuMP.@constraint(wm.model, sum(x_pw[a, :]) == z) - c_3 = JuMP.@constraint(wm.model, lambda[a, 1] <= x_pw[a, 1]) - c_4 = JuMP.@constraint(wm.model, lambda[a, end] <= x_pw[a, end]) - - # Generate a set of uniform flow breakpoints. - breakpoints = range(0.0, stop=JuMP.upper_bound(q), length=pump_breakpoints) - - # Add a constraint for the flow piecewise approximation. - q_lhs = sum(breakpoints[k] * lambda[a, k] for k in 1:pump_breakpoints) - c_5 = JuMP.@constraint(wm.model, q_lhs == q) - - # Add a constraint for the head gain piecewise approximation. - f = _calc_pump_gain_values(collect(breakpoints), curve_fun) - g_lhs = sum(f[k] * lambda[a, k] for k in 1:pump_breakpoints) - c_6 = JuMP.@constraint(wm.model, g_lhs == g) - - # Append the constraint array with the above-generated constraints. - append!(con(wm, n, :head_gain, a), [c_1, c_2, c_3, c_4, c_5, c_6]) - - # Add the adjacency constraints. - for k in 2:pump_breakpoints-1 adjacency = x_pw[a, k-1] + x_pw[a, k] c_7_k = JuMP.@constraint(wm.model, lambda[a, k] <= adjacency) - append!(con(wm, n, :head_gain, a), [c_7_k]) + append!(con(wm, n, :pipe_head_loss, a), [c_7_k]) end end -function constraint_pipe_head_loss_des(wm::LAWaterModel, n::Int, a::Int, alpha::Float64, node_fr::Int, node_to::Int, L::Float64, resistances) - # If the number of breakpoints is not positive, no constraints are added. - pipe_breakpoints = get(wm.ext, :pipe_breakpoints, 0) - if pipe_breakpoints <= 0 return end +"Adds head loss constraint for a design pipe in the `LA` formulation." +function constraint_on_off_pipe_head_loss_des(wm::LAWaterModel, n::Int, a::Int, exponent::Float64, node_fr::Int, node_to::Int, L::Float64, resistances) + # Get the number of breakpoints for the pipe. + n_b = pipe_breakpoints = get(wm.ext, :pipe_breakpoints, 2) # Get common variables and data. - n_b = pipe_breakpoints lambda, x_pw = var(wm, n, :lambda_des_pipe), var(wm, n, :x_pw_des_pipe) h_i, h_j = var(wm, n, :h, node_fr), var(wm, n, :h, node_to) # Initialize flow expression in the head loss relationship. - q_loss_expr = JuMP.AffExpr(0.0) + q_loss_expression = JuMP.AffExpr(0.0) for (r_id, r) in enumerate(resistances) - q, x_res = var(wm, n, :q_des_pipe, a)[r_id], var(wm, n, :x_res, a)[r_id] + q, z = var(wm, n, :q_des_pipe, a)[r_id], var(wm, n, :z_des_pipe, a)[r_id] q_lb, q_ub = JuMP.lower_bound(q), JuMP.upper_bound(q) # Add the first required SOS constraints. lambda_sum = sum(lambda[a, r_id, k] for k in 1:pipe_breakpoints) - c_1 = JuMP.@constraint(wm.model, lambda_sum == x_res) + c_1 = JuMP.@constraint(wm.model, lambda_sum == z) c_2 = JuMP.@constraint(wm.model, lambda[a, r_id, 1] <= x_pw[a, 1]) c_3 = JuMP.@constraint(wm.model, lambda[a, r_id, n_b] <= x_pw[a, n_b-1]) @@ -253,77 +172,35 @@ function constraint_pipe_head_loss_des(wm::LAWaterModel, n::Int, a::Int, alpha:: breakpoints = range(q_lb, stop=q_ub, length=n_b) # Add a constraint for the head loss piecewise approximation. - f = _calc_head_loss_values(collect(breakpoints), alpha) + f = _calc_head_loss_values(collect(breakpoints), exponent) expr_r = r .* sum(f[k] .* lambda[a, r_id, k] for k in 1:pipe_breakpoints) - JuMP.add_to_expression!(q_loss_expr, expr_r) + JuMP.add_to_expression!(q_loss_expression, expr_r) # Add a constraint for the flow piecewise approximation. q_lhs = sum(breakpoints[k] * lambda[a, r_id, k] for k in 1:pipe_breakpoints) c_4 = JuMP.@constraint(wm.model, q_lhs == q) # Append the constraint array with the above-generated constraints. - append!(con(wm, n, :head_loss, a), [c_1, c_2, c_3, c_4]) + append!(con(wm, n, :on_off_pipe_head_loss_des, a), [c_1, c_2, c_3, c_4]) # Add the adjacency constraints. for k in 2:n_b-1 adjacency = x_pw[a, k-1] + x_pw[a, k] c_5_k = JuMP.@constraint(wm.model, lambda[a, r_id, k] <= adjacency) - append!(con(wm, n, :head_loss, a), [c_5_k]) + append!(con(wm, n, :on_off_pipe_head_loss_des, a), [c_5_k]) end end # Add the final SOS and approximation constraints. c_6 = JuMP.@constraint(wm.model, sum(x_pw[a, :]) == 1.0) - c_7 = JuMP.@constraint(wm.model, q_loss_expr == inv(L) * (h_i - h_j)) - append!(con(wm, n, :head_loss, a), [c_6, c_7]) -end - - -function constraint_pipe_head_loss(wm::LAWaterModel, n::Int, a::Int, node_fr::Int, node_to::Int, alpha::Float64, L::Float64, r::Float64) - # If the number of breakpoints is not positive, no constraints are added. - pipe_breakpoints = get(wm.ext, :pipe_breakpoints, 0) - if pipe_breakpoints <= 0 return end - - # Get required variables. - q = var(wm, n, :q_pipe, a) - lambda, x_pw = var(wm, n, :lambda_pipe), var(wm, n, :x_pw_pipe) - h_i, h_j = var(wm, n, :h, node_fr), var(wm, n, :h, node_to) - - # Add the required SOS constraints. - c_1 = JuMP.@constraint(wm.model, sum(lambda[a, :]) == 1.0) - c_2 = JuMP.@constraint(wm.model, sum(x_pw[a, :]) == 1.0) - c_3 = JuMP.@constraint(wm.model, lambda[a, 1] <= x_pw[a, 1]) - c_4 = JuMP.@constraint(wm.model, lambda[a, end] <= x_pw[a, end]) - - # Generate a set of uniform flow breakpoints. - q_lb, q_ub = JuMP.lower_bound(q), JuMP.upper_bound(q) - breakpoints = range(q_lb, stop=q_ub, length=pipe_breakpoints) - - # Add a constraint for the flow piecewise approximation. - q_lhs = sum(breakpoints[k] * lambda[a, k] for k in 1:pipe_breakpoints) - c_5 = JuMP.@constraint(wm.model, q_lhs == q) - - # Add a constraint for the head loss piecewise approximation. - f = _calc_head_loss_values(collect(breakpoints), alpha) - lhs = r * sum(f[k] * lambda[a, k] for k in 1:pipe_breakpoints) - c_6 = JuMP.@constraint(wm.model, lhs == inv(L) * (h_i - h_j)) - - # Append the constraint array with the above-generated constraints. - append!(con(wm, n, :head_loss, a), [c_1, c_2, c_3, c_4, c_5, c_6]) - - # Add the adjacency constraints. - for k in 2:pipe_breakpoints-1 - adjacency = x_pw[a, k-1] + x_pw[a, k] - c_7_k = JuMP.@constraint(wm.model, lambda[a, k] <= adjacency) - append!(con(wm, n, :head_loss, a), [c_7_k]) - end + c_7 = JuMP.@constraint(wm.model, q_loss_expression == inv(L) * (h_i - h_j)) + append!(con(wm, n, :on_off_pipe_head_loss_des, a), [c_6, c_7]) end function objective_owf(wm::LAWaterModel) - # If the number of breakpoints is not positive, no objective is added. - pump_breakpoints = get(wm.ext, :pump_breakpoints, 0) - if pump_breakpoints <= 0 return end + # Get the number of breakpoints for the pump. + pump_breakpoints = get(wm.ext, :pump_breakpoints, 2) # Initialize the objective function. objective = JuMP.AffExpr(0.0) diff --git a/src/form/lrd.jl b/src/form/lrd.jl index 78aba213..c6fd99b9 100644 --- a/src/form/lrd.jl +++ b/src/form/lrd.jl @@ -13,8 +13,8 @@ function _get_head_loss_oa(q::JuMP.VariableRef, q_hat::Float64, alpha::Float64) end -function _get_head_loss_oa_z(q::JuMP.VariableRef, z::Union{JuMP.VariableRef, JuMP.GenericAffExpr}, q_hat::Float64, alpha::Float64) - return q_hat^alpha*z + alpha * q_hat^(alpha - 1.0) * (q - q_hat*z) +function _get_head_loss_oa_binary(q::JuMP.VariableRef, z::Union{JuMP.VariableRef, JuMP.GenericAffExpr}, q_hat::Float64, exponent::Float64) + return q_hat^exponent*z + exponent * q_hat^(exponent - 1.0) * (q - q_hat*z) end @@ -25,157 +25,181 @@ function _get_head_gain_oa(q::JuMP.VariableRef, z::JuMP.VariableRef, q_hat::Floa end -"Pump head gain constraint when the pump status is ambiguous." -function constraint_pump_head_gain(wm::LRDWaterModel, n::Int, a::Int, node_fr::Int, node_to::Int, pc::Array{Float64}) - # If the number of breakpoints is not positive, no constraints are added. - pump_breakpoints = get(wm.ext, :pump_breakpoints, 0) - if pump_breakpoints <= 0 return end +########################################## PIPES ########################################## - # Gather flow, head gain, and convex combination variables. - qp, g, z = var(wm, n, :qp_pump, a), var(wm, n, :g, a), var(wm, n, :z_pump, a) - lambda, x_pw = var(wm, n, :lambda_pump), var(wm, n, :x_pw_pump) - # Add the required SOS constraints. - c_1 = JuMP.@constraint(wm.model, sum(lambda[a, :]) == z) - c_2 = JuMP.@constraint(wm.model, sum(x_pw[a, :]) == z) - c_3 = JuMP.@constraint(wm.model, lambda[a, 1] <= x_pw[a, 1]) - c_4 = JuMP.@constraint(wm.model, lambda[a, end] <= x_pw[a, end]) +function constraint_pipe_head_loss(wm::LRDWaterModel, n::Int, a::Int, node_fr::Int, node_to::Int, exponent::Float64, L::Float64, r::Float64) + # Get the number of breakpoints for the pipe. + num_breakpoints = get(wm.ext, :pipe_breakpoints, 1) - # Add a constraint for the flow piecewise approximation. - breakpoints = range(0.0, stop=JuMP.upper_bound(qp), length=pump_breakpoints) - qp_lhs = sum(breakpoints[k] * lambda[a, k] for k in 1:pump_breakpoints) - c_5 = JuMP.@constraint(wm.model, qp_lhs == qp) + # Get the variable for flow directionality. + y = var(wm, n, :y_pipe, a) - # Append the constraint array. - append!(con(wm, n, :head_gain, a), [c_1, c_2, c_3, c_4, c_5]) + # Get variables for positive flow and head difference. + qp, dhp = var(wm, n, :qp_pipe, a), var(wm, n, :dhp_pipe, a) - for k in 2:pump_breakpoints-1 - # Add the adjacency constraints for piecewise variables. - adjacency = x_pw[a, k-1] + x_pw[a, k] - c_6_k = JuMP.@constraint(wm.model, lambda[a, k] <= adjacency) - append!(con(wm, n, :head_gain, a), [c_6_k]) + # Loop over breakpoints strictly between the lower and upper variable bounds. + for pt in range(0.0, stop = JuMP.upper_bound(qp), length = num_breakpoints+2)[2:end-1] + # Add a linear outer approximation of the convex relaxation at `pt`. + lhs = _get_head_loss_oa_binary(qp, y, pt, exponent) + c = JuMP.@constraint(wm.model, r * lhs <= inv(L) * dhp) + append!(con(wm, n, :pipe_head_loss)[a], [c]) end - # Add head gain outer (i.e., upper) approximations. - for qp_hat in breakpoints - lhs = _get_head_gain_oa(qp, z, qp_hat, pc) - c_7_k = JuMP.@constraint(wm.model, g <= lhs) - append!(con(wm, n, :head_gain, a), [c_7_k]) + # Add linear upper bounds on the above outer approximations. + rhs = r * JuMP.upper_bound(qp)^(exponent - 1.0) * qp + c = JuMP.@constraint(wm.model, inv(L) * dhp <= rhs) + + # Append the :pipe_head_loss constraint array. + append!(con(wm, n, :pipe_head_loss)[a], [c]) + + # Get variables for negative flow and head difference. + qn, dhn = var(wm, n, :qn_pipe, a), var(wm, n, :dhn_pipe, a) + + # Loop over breakpoints strictly between the lower and upper variable bounds. + for pt in range(0.0, stop = JuMP.upper_bound(qn), length = num_breakpoints+2)[2:end-1] + # Add a linear outer approximation of the convex relaxation at `pt`. + lhs = _get_head_loss_oa_binary(qn, 1.0 - y, pt, exponent) + c = JuMP.@constraint(wm.model, r * lhs <= inv(L) * dhn) + append!(con(wm, n, :pipe_head_loss)[a], [c]) end + + # Add linear upper bounds on the above outer approximations. + rhs = r * JuMP.upper_bound(qn)^(exponent - 1.0) * qn + c = JuMP.@constraint(wm.model, inv(L) * dhn <= rhs) + + # Append the :pipe_head_loss constraint array. + append!(con(wm, n, :pipe_head_loss)[a], [c]) end -function constraint_pipe_head_loss_des(wm::LRDWaterModel, n::Int, a::Int, alpha::Float64, node_fr::Int, node_to::Int, L::Float64, resistances) +function constraint_on_off_pipe_head_loss_des(wm::LRDWaterModel, n::Int, a::Int, exponent::Float64, node_fr::Int, node_to::Int, L::Float64, resistances) # If the number of breakpoints is not positive, no constraints are added. - pipe_breakpoints = get(wm.ext, :pipe_breakpoints, 0) - if pipe_breakpoints <= 0 return end + num_breakpoints = get(wm.ext, :pipe_breakpoints, 1) + + # Get design pipe direction and status variable references. + y, z = var(wm, n, :y_des_pipe, a), var(wm, n, :z_des_pipe, a) + + # Get head difference variables for the pipe. + dhp, dhn = var(wm, n, :dhp_des_pipe, a), var(wm, n, :dhn_des_pipe, a) + qp, qn = var(wm, n, :qp_des_pipe, a), var(wm, n, :qn_des_pipe, a) for (r_id, r) in enumerate(resistances) - # Add constraints corresponding to positive outer approximations. - qp, dhp = var(wm, n, :qp_des_pipe, a)[r_id], var(wm, n, :dhp_des_pipe, a) - qp_ub = JuMP.upper_bound(qp) - - for qp_hat in range(0.0, stop=qp_ub, length=pipe_breakpoints) - lhs = r * _get_head_loss_oa(qp, qp_hat, alpha) - c_p = JuMP.@constraint(wm.model, lhs <= inv(L) * dhp) - append!(con(wm, n, :head_loss)[a], [c_p]) + # Get directed flow variables and associated data. + qp_ub, qn_ub = JuMP.upper_bound(qp[r_id]), JuMP.upper_bound(qn[r_id]) + + # Loop over breakpoints strictly between the lower and upper variable bounds. + for pt in range(0.0, stop = qp_ub, length = num_breakpoints+2)[2:end-1] + lhs = r * _get_head_loss_oa_binary(qp[r_id], z[r_id], pt, exponent) + c_1 = JuMP.@constraint(wm.model, lhs <= inv(L) * dhp) + append!(con(wm, n, :on_off_pipe_head_loss_des)[a], [c_1]) end - # Add constraints corresponding to negative outer approximations. - qn, dhn = var(wm, n, :qn_des_pipe, a)[r_id], var(wm, n, :dhn_des_pipe, a) - qn_ub = JuMP.upper_bound(qn) - - for qn_hat in range(0.0, stop=qn_ub, length=pipe_breakpoints) - lhs = r * _get_head_loss_oa(qn, qn_hat, alpha) - c_n = JuMP.@constraint(wm.model, lhs <= inv(L) * dhn) - append!(con(wm, n, :head_loss)[a], [c_n]) + # Loop over breakpoints strictly between the lower and upper variable bounds. + for pt in range(0.0, stop = qn_ub, length = num_breakpoints+2)[2:end-1] + lhs = r * _get_head_loss_oa_binary(qn[r_id], z[r_id], pt, exponent) + c_2 = JuMP.@constraint(wm.model, lhs <= inv(L) * dhn) + append!(con(wm, n, :on_off_pipe_head_loss_des)[a], [c_2]) end end + + # Add linear upper bounds for the positive portion of head loss. + qp_ub = JuMP.upper_bound.(qp) + slopes_p = resistances .* qp_ub.^(exponent - 1.0) + c_3 = JuMP.@constraint(wm.model, inv(L)*dhp <= sum(slopes_p .* qp)) + + # Add linear upper bounds for the negative portion of head loss. + qn_ub = JuMP.upper_bound.(qn) + slopes_n = resistances .* qn_ub.^(exponent - 1.0) + c_4 = JuMP.@constraint(wm.model, inv(L)*dhn <= sum(slopes_n .* qn)) + + # Append the :on_off_pipe_head_loss_des constraint array. + append!(con(wm, n, :on_off_pipe_head_loss_des)[a], [c_3, c_4]) end -function constraint_pipe_head_loss(wm::LRDWaterModel, n::Int, a::Int, node_fr::Int, node_to::Int, alpha::Float64, L::Float64, r::Float64) - # If the number of breakpoints is not positive, no constraints are added. - pipe_breakpoints = get(wm.ext, :pipe_breakpoints, 0) - if pipe_breakpoints <= 0 return end +########################################## PUMPS ########################################## - # Get the variable denoting flow directionality. - y = var(wm, n, :y_pipe, a) - # Add constraints corresponding to positive outer-approximations. - qp, dhp = var(wm, n, :qp_pipe, a), var(wm, n, :dhp_pipe, a) +"Instantiate variables associated with modeling each pump's head gain." +function variable_pump_head_gain(wm::LRDWaterModel; nw::Int=wm.cnw, bounded::Bool=true, report::Bool=true) + # Initialize variables for total hydraulic head gain from a pump. + g = var(wm, nw)[:g_pump] = JuMP.@variable(wm.model, [a in ids(wm, nw, :pump)], + base_name="$(nw)_g_pump", lower_bound=0.0, # Pump gain is nonnegative. + start=comp_start_value(ref(wm, nw, :pump, a), "g_pump_start")) - for qp_hat in range(0.0, stop=JuMP.upper_bound(qp), length=pipe_breakpoints) - lhs = _get_head_loss_oa_z(qp, y, qp_hat, alpha) - c_p = JuMP.@constraint(wm.model, r * lhs <= inv(L) * dhp) - append!(con(wm, n, :head_loss)[a], [c_p]) - end + # Initialize an entry to the solution component dictionary for head gains. + report && sol_component_value(wm, nw, :pump, :g, ids(wm, nw, :pump), g) - # Add constraints corresponding to positive outer-approximations. - qn, dhn = var(wm, n, :qn_pipe, a), var(wm, n, :dhn_pipe, a) + # Get the number of breakpoints for the pump. + pump_breakpoints = get(wm.ext, :pump_breakpoints, 2) - for qn_hat in range(0.0, stop=JuMP.upper_bound(qn), length=pipe_breakpoints) - lhs = _get_head_loss_oa_z(qn, 1.0 - y, qn_hat, alpha) - c_n = JuMP.@constraint(wm.model, r * lhs <= inv(L) * dhn) - append!(con(wm, n, :head_loss)[a], [c_n]) - end + # Create weights involved in convex combination constraints for pumps. + lambda_pump = var(wm, nw)[:lambda_pump] = JuMP.@variable(wm.model, + [a in ids(wm, nw, :pump), k in 1:pump_breakpoints], + base_name="$(nw)_lambda", lower_bound=0.0, upper_bound=1.0, + start=comp_start_value(ref(wm, nw, :pump, a), "lambda_start", k)) + + # Create binary variables involved in convex combination constraints for pumps. + x_pw_pump = var(wm, nw)[:x_pw_pump] = JuMP.@variable(wm.model, + [a in ids(wm, nw, :pump), k in 1:pump_breakpoints-1], base_name="$(nw)_x_pw", + binary=true, start=comp_start_value(ref(wm, nw, :pump, a), "x_pw_start")) end -function constraint_check_valve_head_loss(wm::LRDWaterModel, n::Int, a::Int, node_fr::Int, node_to::Int, L::Float64, r::Float64) - # If the number of breakpoints is not positive, no constraints are added. - pipe_breakpoints = get(wm.ext, :pipe_breakpoints, 0) - if pipe_breakpoints <= 0 return end - - # Get common variables for outer approximation constraints. - y, z = var(wm, n, :y_check_valve, a), var(wm, n, :z_check_valve, a) - qp, dhp = var(wm, n, :qp_check_valve, a), var(wm, n, :dhp_check_valve, a) - - # Add outer approximation constraints. - for qp_hat in range(0.0, stop=JuMP.upper_bound(qp), length=pipe_breakpoints) - lhs_1 = _get_head_loss_oa_z(qp, y, qp_hat, ref(wm, n, :alpha)) - c_p_1 = JuMP.@constraint(wm.model, r * lhs_1 <= inv(L) * dhp) - lhs_2 = _get_head_loss_oa_z(qp, z, qp_hat, ref(wm, n, :alpha)) - c_p_2 = JuMP.@constraint(wm.model, r * lhs_2 <= inv(L) * dhp) - append!(con(wm, n, :head_loss)[a], [c_p_1, c_p_2]) - end -end +"Add constraints associated with modeling a pump's head gain." +function constraint_on_off_pump_head_gain(wm::LRDWaterModel, n::Int, a::Int, node_fr::Int, node_to::Int, pc::Array{Float64}, q_min_active::Float64) + # Get the number of breakpoints for the pump. + pump_breakpoints = get(wm.ext, :pump_breakpoints, 2) + # Gather pump flow, head gain, and status variables. + qp, g, z = var(wm, n, :qp_pump, a), var(wm, n, :g_pump, a), var(wm, n, :z_pump, a) -function constraint_shutoff_valve_head_loss(wm::LRDWaterModel, n::Int, a::Int, node_fr::Int, node_to::Int, L::Float64, r::Float64) - # If the number of breakpoints is not positive, no constraints are added. - pipe_breakpoints = get(wm.ext, :pipe_breakpoints, 0) - if pipe_breakpoints <= 0 return end - - # Get common data for outer approximation constraints. - qp, qn = var(wm, n, :qp_shutoff_valve, a), var(wm, n, :qn_shutoff_valve, a) - dhp, dhn = var(wm, n, :dhp_shutoff_valve, a), var(wm, n, :dhn_shutoff_valve, a) - y, z = var(wm, n, :y_shutoff_valve, a), var(wm, n, :z_shutoff_valve, a) - - # Add outer approximation constraints for positively-directed flow. - for qp_hat in range(0.0, stop=JuMP.upper_bound(qp), length=pipe_breakpoints) - lhs_1 = _get_head_loss_oa_z(qp, y, qp_hat, ref(wm, n, :alpha)) - c_p_1 = JuMP.@constraint(wm.model, L * r * lhs_1 <= dhp) - lhs_2 = _get_head_loss_oa_z(qp, z, qp_hat, ref(wm, n, :alpha)) - c_p_2 = JuMP.@constraint(wm.model, L * r * lhs_2 <= dhp) - append!(con(wm, n, :head_loss)[a], [c_p_1, c_p_2]) + # Gather convex combination variables. + lambda, x_pw = var(wm, n, :lambda_pump), var(wm, n, :x_pw_pump) + + # Add the required SOS constraints. + c_1 = JuMP.@constraint(wm.model, sum(lambda[a, :]) == z) + c_2 = JuMP.@constraint(wm.model, sum(x_pw[a, :]) == z) + c_3 = JuMP.@constraint(wm.model, lambda[a, 1] <= x_pw[a, 1]) + c_4 = JuMP.@constraint(wm.model, lambda[a, end] <= x_pw[a, end]) + + # Add a constraint for the flow piecewise approximation. + breakpoints = range(q_min_active, stop=JuMP.upper_bound(qp), length=pump_breakpoints) + qp_lhs = sum(breakpoints[k] * lambda[a, k] for k in 1:pump_breakpoints) + c_5 = JuMP.@constraint(wm.model, qp_lhs == qp) + + # Add a constraint that lower-bounds the head gain variable. + f = (pc[1] .* breakpoints.^2) .+ (pc[2] .* breakpoints) .+ pc[3] + gain_lb_expr = sum(f[k] .* lambda[a, k] for k in 1:pump_breakpoints) + c_6 = JuMP.@constraint(wm.model, gain_lb_expr <= g) + + # Append the constraint array. + append!(con(wm, n, :on_off_pump_head_gain, a), [c_1, c_2, c_3, c_4, c_5, c_6]) + + for qp_hat in breakpoints + # Add head gain outer (i.e., upper) approximations. + lhs = _get_head_gain_oa(qp, z, qp_hat, pc) + c_7_k = JuMP.@constraint(wm.model, g <= lhs) + append!(con(wm, n, :on_off_pump_head_gain, a), [c_7_k]) end - # Add outer approximation constraints for negatively-directed flow. - for qn_hat in range(0.0, stop=JuMP.upper_bound(qn), length=pipe_breakpoints) - lhs_1 = _get_head_loss_oa_z(qn, 1.0 - y, qn_hat, ref(wm, n, :alpha)) - c_n_1 = JuMP.@constraint(wm.model, L * r * lhs_1 <= dhn) - lhs_2 = _get_head_loss_oa_z(qn, z, qn_hat, ref(wm, n, :alpha)) - c_n_2 = JuMP.@constraint(wm.model, L * r * lhs_2 <= dhn) - append!(con(wm, n, :head_loss)[a], [c_n_1, c_n_2]) + for k in 2:pump_breakpoints-1 + # Add the adjacency constraints for piecewise variables. + adjacency = x_pw[a, k-1] + x_pw[a, k] + c_8_k = JuMP.@constraint(wm.model, lambda[a, k] <= adjacency) + append!(con(wm, n, :on_off_pump_head_gain, a), [c_8_k]) end end +######################################## OBJECTIVES ######################################## + + +"Instantiate the objective associated with the Optimal Water Flow problem." function objective_owf(wm::LRDWaterModel) - # If the number of breakpoints is not positive, no objective is added. - pump_breakpoints = get(wm.ext, :pump_breakpoints, 0) - if pump_breakpoints <= 0 return end + # Get the number of breakpoints for the pump. + pump_breakpoints = get(wm.ext, :pump_breakpoints, 2) # Initialize the objective function. objective = JuMP.AffExpr(0.0) @@ -196,12 +220,11 @@ function objective_owf(wm::LRDWaterModel) head_curve = ref(wm, n, :pump, a)["head_curve"] curve_fun = _get_function_from_head_curve(head_curve) - # Get flow-related variables and data. + # Get pump flow and status variables. qp, z = var(wm, n, :qp_pump, a), var(wm, n, :z_pump, a) - qp_ub = JuMP.upper_bound(qp) # Generate a set of uniform flow and cubic function breakpoints. - breakpoints = range(0.0, stop=qp_ub, length=pump_breakpoints) + breakpoints = range(0.0, stop=JuMP.upper_bound(qp), length=pump_breakpoints) f = _calc_cubic_flow_values(collect(breakpoints), curve_fun) # Get pump efficiency data. diff --git a/src/form/nc.jl b/src/form/nc.jl index f9bf1b65..bb956f2f 100644 --- a/src/form/nc.jl +++ b/src/form/nc.jl @@ -6,93 +6,50 @@ # flow is used in this formulation to define the consumption of power by active pumps). -"Adds head loss constraints for check valves in `NC` formulations." -function constraint_check_valve_head_loss(wm::NCWaterModel, n::Int, a::Int, node_fr::Int, node_to::Int, L::Float64, r::Float64) - # Gather common variables and data. - q, z = var(wm, n, :q_check_valve, a), var(wm, n, :z_check_valve, a) - h_i, h_j = var(wm, n, :h, node_fr), var(wm, n, :h, node_to) - dh_lb = min(0.0, JuMP.lower_bound(h_i) - JuMP.upper_bound(h_j)) - - # There are two possibilities, here: (i) z = 1, in which the check valve is - # open, and the head loss relationship should be met with equality. In this - # case, the constraints below reduce to lhs = 0.0, as would be expected. - # Otherwise, (ii) z = 0, and the head difference must be negative and - # decoupled from the traditional head loss relationship. - lhs = JuMP.@NLexpression(wm.model, inv(L) * (h_i - h_j) - r * head_loss(q)) - c_1 = JuMP.@NLconstraint(wm.model, lhs <= 0.0) - c_2 = JuMP.@NLconstraint(wm.model, lhs >= inv(L) * (1.0 - z) * dh_lb) - - # Append the :head_loss constraint array. - append!(con(wm, n, :head_loss)[a], [c_1, c_2]) -end - - -"Adds head loss constraints for shutoff valves in `NC` formulations." -function constraint_shutoff_valve_head_loss(wm::NCWaterModel, n::Int, a::Int, node_fr::Int, node_to::Int, L::Float64, r::Float64) - # Gather common variables and data. - q, z = var(wm, n, :q_shutoff_valve, a), var(wm, n, :z_shutoff_valve, a) - h_i, h_j = var(wm, n, :h, node_fr), var(wm, n, :h, node_to) - dh_lb = min(0.0, JuMP.lower_bound(h_i) - JuMP.upper_bound(h_j)) - dh_ub = max(0.0, JuMP.upper_bound(h_i) - JuMP.lower_bound(h_j)) - - # There are two possibilities, here: (i) z = 1, in which the shutoff valve - # is open, and the head loss relationship should be met with equality. In - # this case, the constraints below reduce to lhs = 0.0, as would be - # expected. Otherwise, (ii) z = 0, and the head difference must be - # decoupled from the traditional head loss relationship. - lhs = JuMP.@NLexpression(wm.model, (h_i - h_j) - L * r * head_loss(q)) - c_1 = JuMP.@NLconstraint(wm.model, lhs <= (1.0 - z) * dh_ub) - c_2 = JuMP.@NLconstraint(wm.model, lhs >= (1.0 - z) * dh_lb) - - # Append the :head_loss constraint array. - append!(con(wm, n, :head_loss)[a], [c_1, c_2]) -end - - -"Adds head loss constraints for pipes (without check valves) in `NC` formulations." +"Adds head loss constraint for a pipe in the `NC` formulation." function constraint_pipe_head_loss(wm::NCWaterModel, n::Int, a::Int, node_fr::Int, node_to::Int, alpha::Float64, L::Float64, r::Float64) # Gather flow and head variables included in head loss constraints. - q = var(wm, n, :q_pipe, a) - h_i, h_j = var(wm, n, :h, node_fr), var(wm, n, :h, node_to) + 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 = JuMP.@NLconstraint(wm.model, r * head_loss(q) == inv(L) * (h_i - h_j)) - # Append the :head_loss constraint array. - append!(con(wm, n, :head_loss)[a], [c]) + # Append the :pipe_head_loss constraint array. + append!(con(wm, n, :pipe_head_loss)[a], [c]) end -"Add head loss constraints for design pipes (without check valves) in `NC` formulations." -function constraint_pipe_head_loss_des(wm::NCWaterModel, n::Int, a::Int, alpha::Float64, node_fr::Int, node_to::Int, L::Float64, res) +"Adds head loss constraint for a design pipe in the `NC` formulation." +function constraint_on_off_pipe_head_loss_des(wm::NCWaterModel, n::Int, a::Int, exponent::Float64, node_fr::Int, node_to::Int, L::Float64, resistances) # Gather common flow and head variables, as well as design indices. - q, R = var(wm, n, :q_des_pipe, a), 1:length(res) + q, R = var(wm, n, :q_des_pipe, a), 1:length(resistances) h_i, h_j = var(wm, n, :h, node_fr), var(wm, n, :h, node_to) # Add the nonconvex, design-expanded head loss constraint. - lhs = JuMP.@NLexpression(wm.model, sum(res[r] * head_loss(q[r]) for r in R)) + lhs = JuMP.@NLexpression(wm.model, sum(resistances[r] * head_loss(q[r]) for r in R)) c = JuMP.@NLconstraint(wm.model, lhs == inv(L) * (h_i - h_j)) - # Append the :head_loss constraint array. - append!(con(wm, n, :head_loss)[a], [c]) + # Append the :on_off_pipe_head_loss_des constraint array. + append!(con(wm, n, :on_off_pipe_head_loss_des)[a], [c]) end "Adds head gain constraints for pumps in `NC` formulations." -function constraint_pump_head_gain(wm::NCWaterModel, n::Int, a::Int, node_fr::Int, node_to::Int, pc::Array{Float64}) - # Gather common flow and pump variables. - q, g, z = var(wm, n, :q_pump, a), var(wm, n, :g, a), var(wm, n, :z_pump, a) +function constraint_on_off_pump_head_gain(wm::NCWaterModel, n::Int, a::Int, node_fr::Int, node_to::Int, pc::Array{Float64}, q_min_active::Float64) + # Gather pump flow, head gain, and status variables. + q, g, z = var(wm, n, :q_pump, a), var(wm, n, :g_pump, a), var(wm, n, :z_pump, a) # Add constraint equating head gain with respect to the pump curve. c = JuMP.@constraint(wm.model, pc[1]*q^2 + pc[2]*q + pc[3]*z == g) - # Append the :head_gain constraint array. - append!(con(wm, n, :head_gain)[a], [c]) + # Append the :on_off_pump_head_gain constraint array. + append!(con(wm, n, :on_off_pump_head_gain)[a], [c]) end + "Defines the objective for the owf problem is `NC` formulations." function objective_owf(wm::NCWaterModel) - objective = JuMP.AffExpr(0.0) + objective = zero(JuMP.QuadExpr) for (n, nw_ref) in nws(wm) efficiency = 0.85 # TODO: How can the efficiency curve be used? @@ -102,15 +59,11 @@ function objective_owf(wm::NCWaterModel) for (a, pump) in nw_ref[:pump] if haskey(pump, "energy_price") - # Get price data and power-related variables. - price = pump["energy_price"] - q, g = var(wm, n, :q_pump, a), var(wm, n, :g, a) + # Get pump flow and head gain variables. + q, g = var(wm, n, :q_pump, a), var(wm, n, :g_pump, a) # Constrain cost_var and append to the objective expression. - cost_var = JuMP.@variable(wm.model, lower_bound=0.0) - cost_expr = JuMP.@NLexpression(wm.model, coeff * price * g * q) - c = JuMP.@NLconstraint(wm.model, cost_expr <= cost_var) - JuMP.add_to_expression!(objective, cost_var) + JuMP.add_to_expression!(objective, coeff * pump["energy_price"], q, g) else Memento.error(_LOGGER, "No cost given for pump \"$(pump["name"])\"") end diff --git a/src/form/qrd.jl b/src/form/qrd.jl new file mode 100644 index 00000000..d0067395 --- /dev/null +++ b/src/form/qrd.jl @@ -0,0 +1,118 @@ +# Define common QRD (quadratic relaxation- and direction-based) implementations +# of water distribution network constraints, which use directed flow variables. + + +function constraint_pipe_head_loss(wm::QRDWaterModel, n::Int, a::Int, node_fr::Int, node_to::Int, exponent::Float64, L::Float64, r::Float64) + # Get the number of breakpoints for the pipe. + num_breakpoints = get(wm.ext, :pipe_breakpoints, 1) + + # 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) + + # Loop over breakpoints strictly between the lower and upper variable bounds. + for pt in range(0.0, stop = JuMP.upper_bound(qp), length = num_breakpoints+2)[2:end-1] + # Add a linear outer approximation of the convex relaxation at `pt`. + lhs = _get_head_loss_oa_binary(qp, y, pt, exponent) + c = JuMP.@constraint(wm.model, r * lhs <= inv(L) * dhp) + append!(con(wm, n, :pipe_head_loss)[a], [c]) + end + + # Add linear upper bounds on the above outer approximations. + rhs = r * JuMP.upper_bound(qp)^(exponent - 1.0) * qp + c = JuMP.@constraint(wm.model, inv(L) * dhp <= rhs) + + # Append the :pipe_head_loss constraint array. + append!(con(wm, n, :pipe_head_loss)[a], [c]) + + # Get variables for negative flow and head difference. + qn, dhn = var(wm, n, :qn_pipe, a), var(wm, n, :dhn_pipe, a) + + # Loop over breakpoints strictly between the lower and upper variable bounds. + for pt in range(0.0, stop = JuMP.upper_bound(qn), length = num_breakpoints+2)[2:end-1] + # Add a linear outer approximation of the convex relaxation at `pt`. + lhs = _get_head_loss_oa_binary(qn, 1.0 - y, pt, exponent) + c = JuMP.@constraint(wm.model, r * lhs <= inv(L) * dhn) + append!(con(wm, n, :pipe_head_loss)[a], [c]) + end + + # Add linear upper bounds on the above outer approximations. + rhs = r * JuMP.upper_bound(qn)^(exponent - 1.0) * qn + c = JuMP.@constraint(wm.model, inv(L) * dhn <= rhs) + + # Append the :pipe_head_loss constraint array. + append!(con(wm, n, :pipe_head_loss)[a], [c]) +end + + +"Pump head gain constraint when the pump status is ambiguous." +function constraint_on_off_pump_head_gain(wm::QRDWaterModel, n::Int, a::Int, node_fr::Int, node_to::Int, pc::Array{Float64}, q_min_active::Float64) + # Gather pump flow, head gain, and status variables. + qp, g, z = var(wm, n, :qp_pump, a), var(wm, n, :g_pump, a), var(wm, n, :z_pump, a) + + # Define the (relaxed) head gain relationship for the pump. + c = JuMP.@constraint(wm.model, g == pc[1]*qp^2 + pc[2]*qp + pc[3]*z) + append!(con(wm, n, :on_off_pump_head_gain, a), [c]) +end + + +function objective_wf(wm::QRDWaterModel) + JuMP.set_objective_sense(wm.model, _MOI.FEASIBILITY_SENSE) +end + + +function objective_owf(wm::QRDWaterModel) + # Initialize the objective function. + objective = zero(JuMP.QuadExpr) + + for (n, nw_ref) in nws(wm) + # Get common constant parameters. + rho = 1000.0 # Water density (kilogram per cubic meter). + gravity = 9.80665 # Gravitational acceleration (meter per second squared). + constant = rho * gravity * ref(wm, n, :time_step) + + for (a, pump) in nw_ref[:pump] + q_min_active = get(pump, "q_min_active", _q_eps) + + if haskey(pump, "energy_price") + # Get price and pump curve data. + price = pump["energy_price"] + head_curve = ref(wm, n, :pump, a)["head_curve"] + curve_fun = _get_function_from_head_curve(head_curve) + + # Get flow-related variables and data. + z = var(wm, n, :z_pump, a) + qp, g = var(wm, n, :qp_pump, a), var(wm, n, :g_pump, a) + points = collect(range(q_min_active, stop=JuMP.upper_bound(qp), length=50)) + + # Get pump efficiency data. + if haskey(pump, "efficiency_curve") + # Use the efficiency at the midpoint of the pump flow bounds. + eff_curve = pump["efficiency_curve"] + eff = _calc_efficiencies(points, eff_curve) + else + eff = pump["efficiency"] + end + + # Compute discrete costs from existing points. + flows_cubed = _calc_cubic_flow_values(points, curve_fun) + costs = (constant*price) .* inv.(eff) .* flows_cubed + + # Fit a quadratic function to the above discrete costs. + LsqFit.@. func(x, p) = p[1]*x*x + p[2]*x + p[3] + fit = LsqFit.curve_fit(func, points, costs, zeros(length(costs))) + coeffs = LsqFit.coef(fit) + + # Add the cost corresponding to the current pump's operation. + JuMP.add_to_expression!(objective, coeffs[1], qp, qp) + JuMP.add_to_expression!(objective, coeffs[2] * qp + coeffs[3] * z) + else + Memento.error(_LOGGER, "No cost given for pump \"$(pump["name"])\"") + end + end + end + + return JuMP.@objective(wm.model, _MOI.MIN_SENSE, objective) +end diff --git a/src/form/undirected_flow.jl b/src/form/undirected_flow.jl index de383303..83666bc3 100644 --- a/src/form/undirected_flow.jl +++ b/src/form/undirected_flow.jl @@ -3,10 +3,11 @@ # When q is nonnegative, flow is assumed to travel from i to j. When q is # negative, flow is assumed to travel from j to i. + "Create flow-related variables common to all directed flow models for edge-type components." function variable_flow(wm::AbstractUndirectedModel; nw::Int=wm.cnw, bounded::Bool=true, report::Bool=true) - for name in ["check_valve", "pipe", "pressure_reducing_valve", "pump", "shutoff_valve"] - # Create directed flow (`qp` and `qn`) variables for each component. + for name in ["pipe", "pump", "regulator", "short_pipe", "valve"] + # Create undirected flow variables for each component. variable_component_flow(wm, name; nw=nw, bounded=bounded, report=report) end @@ -70,155 +71,205 @@ function variable_flow_des_common(wm::AbstractUndirectedModel; nw::Int=wm.cnw, b # Initialize the solution reporting data structures. report && sol_component_value(wm, nw, :des_pipe, :q, ids(wm, nw, :des_pipe), q) +end - # Create resistance binary variables. - variable_resistance(wm, nw=nw) + +function constraint_pipe_flow(wm::AbstractUndirectedModel, n::Int, a::Int) + # By default, there are no constraints, here. +end + + +function constraint_pipe_head(wm::AbstractUndirectedModel, n::Int, a::Int, node_fr::Int, node_to::Int) + # Get head variables for from and to nodes. + h_i, h_j = var(wm, n, :h, node_fr), var(wm, n, :h, node_to) + + # For pipes, the differences must satisfy lower and upper bounds. + dh_lb = JuMP.lower_bound(h_i) - JuMP.upper_bound(h_j) + c_1 = JuMP.@constraint(wm.model, h_i - h_j >= dh_lb) + dh_ub = JuMP.upper_bound(h_i) - JuMP.lower_bound(h_j) + c_2 = JuMP.@constraint(wm.model, h_i - h_j <= dh_ub) + + # Append the constraint array. + append!(con(wm, n, :pipe_head, a), [c_1, c_2]) end "Constrain flow variables, based on design selections, in undirected flow formulations." -function constraint_resistance_selection_des(wm::AbstractUndirectedModel, n::Int, a::Int, pipe_resistances) - c = JuMP.@constraint(wm.model, sum(var(wm, n, :x_res, a)) == 1.0) - append!(con(wm, n, :head_loss)[a], [c]) +function constraint_on_off_pipe_flow_des(wm::AbstractUndirectedModel, n::Int, a::Int, resistances) + # Get design pipe status variable references. + z = var(wm, n, :z_des_pipe, a) - for r in 1:length(pipe_resistances) - q_des_pipe = var(wm, n, :q_des_pipe, a)[r] - x_res = var(wm, n, :x_res, a)[r] + # Ensure that only one flow can be nonnegative per solution. + c_1 = JuMP.@constraint(wm.model, sum(z) == 1.0) + append!(con(wm, n, :on_off_pipe_flow_des)[a], [c_1]) - q_des_pipe_lb = JuMP.lower_bound(q_des_pipe) - c_lb = JuMP.@constraint(wm.model, q_des_pipe >= q_des_pipe_lb * x_res) + for r_id in 1:length(resistances) + # Get directed flow variables and associated data. + q = var(wm, n, :q_des_pipe, a)[r_id] + q_lb, q_ub = JuMP.lower_bound(q), JuMP.upper_bound(q) - q_des_pipe_ub = JuMP.upper_bound(q_des_pipe) - c_ub = JuMP.@constraint(wm.model, q_des_pipe <= q_des_pipe_ub * x_res) + # Constraint the pipes based on direction and construction status. + c_2 = JuMP.@constraint(wm.model, q >= q_lb * z[r_id]) + c_3 = JuMP.@constraint(wm.model, q <= q_ub * z[r_id]) - append!(con(wm, n, :head_loss)[a], [c_lb, c_ub]) + # Append the :on_off_pipe_flow_des constraint array. + append!(con(wm, n, :on_off_pipe_flow_des)[a], [c_2, c_3]) end end -function constraint_check_valve_common(wm::AbstractUndirectedModel, n::Int, a::Int, node_fr::Int, node_to::Int) - # Get flow and check valve status variables. - q, z = var(wm, n, :q_check_valve, a), var(wm, n, :z_check_valve, a) +"Constrain head variables in undirected flow formulations." +function constraint_on_off_pipe_head_des(wm::AbstractUndirectedModel, n::Int, a::Int, node_fr::Int, node_to::Int) + # By default, there are no constraints, here. +end - # If the check valve is closed, flow must be zero. - c_1 = JuMP.@constraint(wm.model, q <= JuMP.upper_bound(q) * z) - # Get head variables for from and to nodes. +function constraint_on_off_regulator_flow(wm::AbstractUndirectedModel, n::Int, a::Int, q_min_active::Float64) + # Get flow and regulator status variables. + q, z = var(wm, n, :q_regulator, a), var(wm, n, :z_regulator, a) + + # If the regulator is closed, flow must be zero. + q_lb, q_ub = max(JuMP.lower_bound(q), q_min_active), JuMP.upper_bound(q) + c_1 = JuMP.@constraint(wm.model, q >= q_lb * z) + c_2 = JuMP.@constraint(wm.model, q <= q_ub * z) + + # Append the constraint array. + append!(con(wm, n, :on_off_regulator_flow, a), [c_1, c_2]) +end + + +function constraint_on_off_regulator_head(wm::AbstractUndirectedModel, n::Int, a::Int, node_fr::Int, node_to::Int, head_setting::Float64) + # Get regulator status variable. + z = var(wm, n, :z_regulator, a) + + # Get head variables for from and to nodes (i.e., `i` and `j`). h_i, h_j = var(wm, n, :h, node_fr), var(wm, n, :h, node_to) - # When the check valve is open, negative head loss is not possible. - dh_lb = JuMP.lower_bound(h_i) - JuMP.upper_bound(h_j) - c_2 = JuMP.@constraint(wm.model, h_i - h_j >= (1.0 - z) * min(0.0, dh_lb)) + # When the pressure reducing valve is open, the head at node j is predefined. + h_lb, h_ub = JuMP.lower_bound(h_j), JuMP.upper_bound(h_j) + c_1 = JuMP.@constraint(wm.model, h_j >= (1.0 - z) * h_lb + z * head_setting) + c_2 = JuMP.@constraint(wm.model, h_j <= (1.0 - z) * h_ub + z * head_setting) - # When the check valve is closed, positive head loss is not possible. + # When the pressure reducing valve is open, the head loss is nonnegative. + dh_lb = JuMP.lower_bound(h_i) - JuMP.upper_bound(h_j) + c_3 = JuMP.@constraint(wm.model, h_i - h_j >= dh_lb * (1.0 - z)) dh_ub = JuMP.upper_bound(h_i) - JuMP.lower_bound(h_j) - c_3 = JuMP.@constraint(wm.model, h_i - h_j <= z * dh_ub) + c_4 = JuMP.@constraint(wm.model, h_i - h_j <= dh_ub * z) # Append the constraint array. - append!(con(wm, n, :check_valve, a), [c_1, c_2, c_3]) + append!(con(wm, n, :on_off_regulator_head, a), [c_1, c_2, c_3, c_4]) end -function constraint_sv_common(wm::AbstractUndirectedModel, n::Int, a::Int, node_fr::Int, node_to::Int) - # Get flow and shutoff valve status variables. - q, z = var(wm, n, :q_shutoff_valve, a), var(wm, n, :z_shutoff_valve, a) +function constraint_on_off_valve_flow(wm::AbstractUndirectedModel, n::Int, a::Int) + # Get flow and valve status variables. + q, z = var(wm, n, :q_valve, a), var(wm, n, :z_valve, a) - # If the shutoff valve is open, flow must be appreciably nonnegative. - c_1 = JuMP.@constraint(wm.model, q <= JuMP.upper_bound(q) * z) - c_2 = JuMP.@constraint(wm.model, q >= JuMP.lower_bound(q) * z) + # If the valve is closed, flow must be zero. + q_lb, q_ub = JuMP.lower_bound(q), JuMP.upper_bound(q) + c_1 = JuMP.@constraint(wm.model, q >= q_lb * z) + c_2 = JuMP.@constraint(wm.model, q <= q_ub * z) # Append the constraint array. - append!(con(wm, n, :sv, a), [c_1, c_2]) + append!(con(wm, n, :on_off_valve_flow, a), [c_1, c_2]) end -function constraint_prv_common(wm::AbstractUndirectedModel, n::Int, a::Int, node_fr::Int, node_to::Int, h_prv::Float64) - # Get flow and pressure reducing valve status variables. - q = var(wm, n, :q_pressure_reducing_valve, a) - z = var(wm, n, :z_pressure_reducing_valve, a) - - # If the pressure reducing valve is open, flow must be appreciably nonnegative. - c_1 = JuMP.@constraint(wm.model, q <= JuMP.upper_bound(q) * z) - c_2 = JuMP.@constraint(wm.model, q >= _q_eps * z) +function constraint_on_off_valve_head(wm::AbstractUndirectedModel, n::Int, a::Int, node_fr::Int, node_to::Int) + # Get flow and valve status variables. + q, z = var(wm, n, :q_valve, a), var(wm, n, :z_valve, a) # Get head variables for from and to nodes. h_i, h_j = var(wm, n, :h, node_fr), var(wm, n, :h, node_to) - # When the pressure reducing valve is open, the head at node j is predefined. - h_lb, h_ub = JuMP.lower_bound(h_j), JuMP.upper_bound(h_j) - c_3 = JuMP.@constraint(wm.model, h_j >= (1.0 - z) * h_lb + z * h_prv) - c_4 = JuMP.@constraint(wm.model, h_j <= (1.0 - z) * h_ub + z * h_prv) - - # When the pressure reducing valve is open, the head loss is nonnegative. + # When the valve is open, negative head loss is not possible. dh_lb = JuMP.lower_bound(h_i) - JuMP.upper_bound(h_j) + c_1 = JuMP.@constraint(wm.model, h_i - h_j >= (1.0 - z) * dh_lb) + + # When the valve is closed, positive head loss is not possible. dh_ub = JuMP.upper_bound(h_i) - JuMP.lower_bound(h_j) - c_5 = JuMP.@constraint(wm.model, h_i - h_j <= dh_ub * z) - c_6 = JuMP.@constraint(wm.model, h_i - h_j >= dh_lb * (1.0 - z)) + c_2 = JuMP.@constraint(wm.model, h_i - h_j <= (1.0 - z) * dh_ub) + + # Append the constraint array. + append!(con(wm, n, :on_off_valve_head, a), [c_1, c_2]) +end + + +function constraint_on_off_pump_flow(wm::AbstractUndirectedModel, n::Int, a::Int, q_min_active::Float64) + # Get pump status variable. + q, z = var(wm, n, :q_pump, a), var(wm, n, :z_pump, a) + + # If the pump is inactive, flow must be zero. + q_ub = JuMP.upper_bound(q) + c_1 = JuMP.@constraint(wm.model, q >= q_min_active * z) + c_2 = JuMP.@constraint(wm.model, q <= q_ub * z) # Append the constraint array. - append!(con(wm, n, :prv, a), [c_1, c_2, c_3, c_4, c_5, c_6]) + append!(con(wm, n, :on_off_pump_flow, a), [c_1, c_2]) end -function constraint_pump_common(wm::AbstractUndirectedModel, n::Int, a::Int, node_fr::Int, node_to::Int, pc::Array{Float64}) - # Gather common variables. - q, g, z = var(wm, n, :q_pump, a), var(wm, n, :g, a), var(wm, n, :z_pump, a) +function constraint_on_off_pump_head(wm::AbstractUndirectedModel, n::Int, a::Int, node_fr::Int, node_to::Int) + # Get head variables for from and to nodes. h_i, h_j = var(wm, n, :h, node_fr), var(wm, n, :h, node_to) - # If the pump is off, the flow along the pump must be zero. - c_1 = JuMP.@constraint(wm.model, q <= JuMP.upper_bound(q) * z) - c_2 = JuMP.@constraint(wm.model, q >= _q_eps * z) + # Get pump status variable. + g, z = var(wm, n, :g_pump, a), var(wm, n, :z_pump, a) - # If the pump is off, decouple the head difference relationship. - dhn_lb = min(0.0, JuMP.lower_bound(h_j) - JuMP.upper_bound(h_i)) - c_3 = JuMP.@constraint(wm.model, h_j - h_i - g >= dhn_lb * (1.0 - z)) - dhn_ub = max(0.0, JuMP.upper_bound(h_j) - JuMP.lower_bound(h_i)) - c_4 = JuMP.@constraint(wm.model, h_j - h_i - g <= dhn_ub * (1.0 - z)) + # If the pump is off, decouple the head difference relationship. If the pump is on, + # ensure the head difference is equal to the pump's head gain (i.e., `g`). + dhn_ub = JuMP.upper_bound(h_j) - JuMP.lower_bound(h_i) + dhn_lb = JuMP.lower_bound(h_j) - JuMP.upper_bound(h_i) + c_1 = JuMP.@constraint(wm.model, h_j - h_i <= g + dhn_ub * (1.0 - z)) + c_2 = JuMP.@constraint(wm.model, h_j - h_i >= g + dhn_lb * (1.0 - z)) # Append the constraint array. - append!(con(wm, n, :pump, a), [c_1, c_2, c_3, c_4]) + append!(con(wm, n, :on_off_pump_head, a), [c_1, c_2]) end -function constraint_pipe_common(wm::AbstractUndirectedModel, n::Int, a::Int, node_fr::Int, node_to::Int, alpha::Float64, L::Float64, r::Float64) - # For undirected formulations, there are no constraints, here. + +function constraint_short_pipe_flow(wm::AbstractUndirectedModel, n::Int, a::Int) + # By default, there are no constraints, here. +end + + +function constraint_short_pipe_head(wm::AbstractUndirectedModel, n::Int, a::Int, node_fr::Int, node_to::Int) + # Get head variables for from and to nodes. + h_i, h_j = var(wm, n, :h, node_fr), var(wm, n, :h, node_to) + + # For short pipes, the heads at adjacent nodes are equal. + c = JuMP.@constraint(wm.model, h_i - h_j == 0.0) + + # Append the constraint array. + append!(con(wm, n, :short_pipe_head, a), [c]) end -function constraint_node_directionality( - wm::AbstractUndirectedModel, n::Int, i::Int, check_valve_fr::Array{Int64,1}, - check_valve_to::Array{Int64,1}, pipe_fr::Array{Int64,1}, pipe_to::Array{Int64,1}, - des_pipe_fr::Array{Int64,1}, des_pipe_to::Array{Int64,1}, pump_fr::Array{Int64,1}, - pump_to::Array{Int64,1}, pressure_reducing_valve_fr::Array{Int64,1}, - pressure_reducing_valve_to::Array{Int64,1}, shutoff_valve_fr::Array{Int64,1}, - shutoff_valve_to::Array{Int64,1}) +function constraint_intermediate_directionality( + wm::AbstractUndirectedModel, n::Int, i::Int, pipe_fr::Array{Int64,1}, + pipe_to::Array{Int64,1}, des_pipe_fr::Array{Int64,1}, des_pipe_to::Array{Int64,1}, + pump_fr::Array{Int64,1}, pump_to::Array{Int64,1}, regulator_fr::Array{Int64,1}, + regulator_to::Array{Int64,1}, short_pipe_fr::Array{Int64,1}, + short_pipe_to::Array{Int64,1}, valve_fr::Array{Int64,1}, valve_to::Array{Int64,1}) # For undirected formulations, there are no constraints, here. end function constraint_sink_directionality( - wm::AbstractUndirectedModel, n::Int, i::Int, check_valve_fr::Array{Int64,1}, - check_valve_to::Array{Int64,1}, pipe_fr::Array{Int64,1}, pipe_to::Array{Int64,1}, - des_pipe_fr::Array{Int64,1}, des_pipe_to::Array{Int64,1}, pump_fr::Array{Int64,1}, - pump_to::Array{Int64,1}, pressure_reducing_valve_fr::Array{Int64,1}, - pressure_reducing_valve_to::Array{Int64,1}, shutoff_valve_fr::Array{Int64,1}, - shutoff_valve_to::Array{Int64,1}) + wm::AbstractUndirectedModel, n::Int, i::Int, pipe_fr::Array{Int64,1}, + pipe_to::Array{Int64,1}, des_pipe_fr::Array{Int64,1}, des_pipe_to::Array{Int64,1}, + pump_fr::Array{Int64,1}, pump_to::Array{Int64,1}, regulator_fr::Array{Int64,1}, + regulator_to::Array{Int64,1}, short_pipe_fr::Array{Int64,1}, + short_pipe_to::Array{Int64,1}, valve_fr::Array{Int64,1}, valve_to::Array{Int64,1}) # For undirected formulations, there are no constraints, here. end function constraint_source_directionality( - wm::AbstractUndirectedModel, n::Int, i::Int, check_valve_fr::Array{Int64,1}, - check_valve_to::Array{Int64,1}, pipe_fr::Array{Int64,1}, pipe_to::Array{Int64,1}, - des_pipe_fr::Array{Int64,1}, des_pipe_to::Array{Int64,1}, pump_fr::Array{Int64,1}, - pump_to::Array{Int64,1}, pressure_reducing_valve_fr::Array{Int64,1}, - pressure_reducing_valve_to::Array{Int64,1}, shutoff_valve_fr::Array{Int64,1}, - shutoff_valve_to::Array{Int64,1}) + wm::AbstractUndirectedModel, n::Int, i::Int, pipe_fr::Array{Int64,1}, + pipe_to::Array{Int64,1}, des_pipe_fr::Array{Int64,1}, des_pipe_to::Array{Int64,1}, + pump_fr::Array{Int64,1}, pump_to::Array{Int64,1}, regulator_fr::Array{Int64,1}, + regulator_to::Array{Int64,1}, short_pipe_fr::Array{Int64,1}, + short_pipe_to::Array{Int64,1}, valve_fr::Array{Int64,1}, valve_to::Array{Int64,1}) # For undirected formulations, there are no constraints, here. end - -function constraint_flow_direction_selection_des(wm::AbstractUndirectedModel, n::Int, a::Int, pipe_resistances) end -function constraint_head_loss_ub_cv(wm::AbstractUndirectedModel, n::Int, a::Int, alpha::Float64, L::Float64, r::Float64) end -function constraint_shutoff_valve_head_loss_ub(wm::AbstractUndirectedModel, n::Int, a::Int, alpha::Float64, L::Float64, r::Float64) end -function constraint_pipe_head_loss_ub_des(wm::AbstractUndirectedModel, n::Int, a::Int, alpha, len, pipe_resistances) end -function constraint_pipe_head_loss_ub(wm::AbstractUndirectedModel, n::Int, a::Int, alpha, len, r_max) end -function constraint_pump_head_gain_lb(wm::AbstractUndirectedModel, n::Int, a::Int, node_fr::Int, node_to::Int, pc::Array{Float64}) end diff --git a/src/io/common.jl b/src/io/common.jl index bb2f068c..efe71246 100644 --- a/src/io/common.jl +++ b/src/io/common.jl @@ -15,6 +15,9 @@ function parse_file(path::String) error("\"$(path)\" is not a valid file type.") end + # Perform data corrections. + _correct_flow_bounds!(network_data) + return network_data end diff --git a/src/io/epanet.jl b/src/io/epanet.jl index 235c10db..399baa0e 100644 --- a/src/io/epanet.jl +++ b/src/io/epanet.jl @@ -105,7 +105,7 @@ function parse_epanet(filename::String) _read_pattern!(data) # Parse [JUNCTIONS] section. - _read_junction!(data) + _read_demand!(data) # Parse [RESERVOIRS] section. _read_reservoir!(data) @@ -119,8 +119,8 @@ function parse_epanet(filename::String) # Add a data structure mapping names to node indices. _add_node_map!(data) - # Update time series information for junctions. - _update_junction_ts!(data) + # Update time series information for demands. + _update_demand_ts!(data) # Update time series information for reservoirs. _update_reservoir_ts!(data) @@ -128,11 +128,17 @@ function parse_epanet(filename::String) # Parse [PIPES] section. _read_pipe!(data) + # Set up a dictionary containing valve objects. + data["valve"] = Dict{String,Any}() + + # Set up a dictionary containing short pipe objects. + data["short_pipe"] = Dict{String,Any}() + # Parse [PUMPS] section. _read_pump!(data) # Parse [VALVES] section. - _read_valve!(data) + _read_regulator!(data) # Parse [ENERGY] section. _read_energy!(data) @@ -168,56 +174,6 @@ function parse_epanet(filename::String) end -""" - epanet_to_watermodels!(epanet_data; import_all=false) - -Converts data parsed from an EPANET file, passed by `epanet_data` into a format suitable for -internal WaterModels use. Imports all data from the EPANET file if `import_all` is true. -""" -function epanet_to_watermodels!(data::Dict{String,<:Any}; import_all::Bool = false) - edge_index, node_index = 0, maximum([x["index"] for (i, x) in data["node"]]) + 1 - - # Determine the starting index of new edges to be added in the network. - for type in ["pipe", "pressure_reducing_valve", "pump"] - if length(data[type]) > 0 - max_component_id = maximum([x["index"] for (i, x) in data[type]]) - edge_index = max(edge_index, max_component_id + 1) - end - end - - # Modify the network for standard modeling of tanks. - for (i, tank) in data["tank"] - # Create a new node, which will be connected to the tank with a shutoff valve. - old_index = string(tank["node"]) - node = deepcopy(data["node"][old_index]) - node["index"], node["name"] = node_index, string(node_index) - data["node"][string(node_index)] = node - # adjust the node source-id for the intermediate node - data["node"][old_index]["source_id"] = AbstractString["node", "$(old_index)"] - - # Instantiate the properties that define the auxiliary pipe. - pipe = Dict{String,Any}("name" => string(edge_index), "status" => 1) - pipe["source_id"] = ["pipe", string(edge_index)] - pipe["node_fr"], pipe["node_to"] = tank["node"], node_index - pipe["length"], pipe["diameter"], pipe["flow_direction"] = 0.0, 1.0, UNKNOWN - pipe["has_check_valve"], pipe["has_shutoff_valve"] = false, true - pipe["minor_loss"], pipe["roughness"] = 0.0, 100.0 - data["pipe"][string(edge_index)] = pipe - - # Set the tank node index to the index of the dummy node. - tank["node"] = node_index - - # Update the auxiliary node and edge indices. - node_index, edge_index = node_index + 1, edge_index + 1 - end - - for (i, junction) in data["junction"] - if isapprox(junction["demand"], 0.0, atol=1.0e-7) - delete!(data["junction"], i) - haskey(data, "time_series") && delete!(data["time_series"]["junction"], i) - end - end -end function _update_time_series!(data::Dict{String,<:Any}) @@ -226,7 +182,7 @@ function _update_time_series!(data::Dict{String,<:Any}) data["time_step"] = convert(Float64, time_step) if num_steps >= 1 && keys(data["pattern"]) != ["1"] - for type in ["junction", "node", "pump"] + for type in ["demand", "node", "pump"] length(data["time_series"][type]) == 0 && delete!(data["time_series"], type) end @@ -261,18 +217,18 @@ function _transform_component_indices(components::Dict{String,<:Any}) end -function _update_junction_ts!(data::Dict{String,<:Any}) - # Create a temporary dictionary representing the junction time series. - junction_ts = Dict{String,Any}() +function _update_demand_ts!(data::Dict{String,<:Any}) + # Create a temporary dictionary representing the demand time series. + demand_ts = Dict{String,Any}() - # Ensure that junction time series data use new junction indices. - for (i, junction) in data["time_series"]["junction"] - key = findfirst(x -> x["source_id"][2] == i, data["junction"]) - junction_ts[key] = junction + # Ensure that demand time series data use new demand indices. + for (i, demand) in data["time_series"]["demand"] + key = findfirst(x -> x["source_id"][2] == i, data["demand"]) + demand_ts[key] = demand end - # Update the junction entry in the time series dictionary. - data["time_series"]["junction"] = junction_ts + # Update the demand entry in the time series dictionary. + data["time_series"]["demand"] = demand_ts end @@ -295,7 +251,7 @@ end function _add_nodes!(data::Dict{String,<:Any}) # Define EPANET nodal types that should be attached to nodes. - comp_types = ["junction", "reservoir", "tank"] + comp_types = ["demand", "reservoir", "tank"] # Obtain the original EPANET indices for all nodal components. comp_names = vcat([collect(keys(data[t])) for t in comp_types]...) @@ -654,11 +610,11 @@ function _update_pump_energy!(data::Dict{String,<:Any}) end -"Parse and store junction data from an EPANET-formatted data dictionary." -function _read_junction!(data::Dict{String,<:Any}) - # Initialize dictionaries associated with junctions. - data["junction"] = Dict{String,Dict{String,Any}}() - data["time_series"]["junction"] = Dict{String,Any}() +"Parse and store demand data from an EPANET-formatted data dictionary." +function _read_demand!(data::Dict{String,<:Any}) + # Initialize dictionaries associated with demands. + data["demand"] = Dict{String,Dict{String,Any}}() + data["time_series"]["demand"] = Dict{String,Any}() # Initialize a temporary index to be updated while parsing. index::Int64 = 0 @@ -668,23 +624,23 @@ function _read_junction!(data::Dict{String,<:Any}) current = split(split(line, ";")[1]) length(current) == 0 && continue # Skip. - # Initialize the junction entry to be added. - junction = Dict{String,Any}("name" => current[1], "status" => 1) - junction["source_id"] = ["junction", current[1]] - junction["dispatchable"] = false + # Initialize the demand entry to be added. + demand = Dict{String,Any}("name" => current[1], "status" => 1) + demand["source_id"] = ["demand", current[1]] + demand["dispatchable"] = false - # Parse the elevation of the junction (in meters). + # Parse the elevation of the demand (in meters). if data["flow_units"] == "LPS" || data["flow_units"] == "CMH" # Retain the original value (in meters). - junction["elevation"] = parse(Float64, current[2]) + demand["elevation"] = parse(Float64, current[2]) elseif data["flow_units"] == "GPM" # If gallons per minute... # Convert elevation from feet to meters. - junction["elevation"] = 0.3048 * parse(Float64, current[2]) + demand["elevation"] = 0.3048 * parse(Float64, current[2]) else Memento.error(_LOGGER, "Could not find a valid \"units\" option type.") end - # Parse the demand at the junction (in cubic meters per second). + # Parse the demand at the demand (in cubic meters per second). if length(current) > 2 # A value for demand must exist. # Calculate the unscaled (unit unknown) demand. demand_unscaled = parse(Float64, current[3]) * data["demand_multiplier"] @@ -692,40 +648,40 @@ function _read_junction!(data::Dict{String,<:Any}) # Convert the unscaled demand to cubic meters per second. if data["flow_units"] == "LPS" # If liters per second... # Convert from liters per second to cubic meters per second. - junction["demand"] = 1.0e-3 * demand_unscaled + demand["flow_rate"] = 1.0e-3 * demand_unscaled elseif data["flow_units"] == "CMH" # If cubic meters per hour... # Convert from cubic meters per hour to cubic meters per second. - junction["demand"] = inv(3600.0) * demand_unscaled + demand["flow_rate"] = inv(3600.0) * demand_unscaled elseif data["flow_units"] == "GPM" # If gallons per minute... # Convert from gallons per minute to cubic meters per second. - junction["demand"] = 6.30902e-5 * demand_unscaled + demand["flow_rate"] = 6.30902e-5 * demand_unscaled end else # No value for demand exists in the line. - junction["status"] = 0 # Set the status of the junction to zero. - junction["demand"] = 0.0 # Initialize demand to zero. + demand["status"] = 0 # Set the status of the demand to zero. + demand["flow_rate"] = 0.0 # Initialize demand to zero. end - # Parse the name of the pattern used to scale demand at the junction. + # Parse the name of the pattern used to scale demand at the demand. pattern = length(current) > 3 ? current[4] : data["default_pattern"] # Scale the parsed demand by the specified pattern, if it exists. if pattern != nothing && length(data["pattern"][pattern]) > 1 - demand = junction["demand"] .* data["pattern"][pattern] - entry = Dict{String,Array{Float64}}("demand" => demand) - data["time_series"]["junction"][current[1]] = entry + flow_rate = demand["flow_rate"] .* data["pattern"][pattern] + entry = Dict{String,Array{Float64}}("flow_rate" => flow_rate) + data["time_series"]["demand"][current[1]] = entry elseif pattern != nothing && pattern != "1" - junction["demand"] *= data["pattern"][pattern][1] + demand["flow_rate"] *= data["pattern"][pattern][1] end # Add a temporary index to be used in the data dictionary. - junction["index"] = string(index += 1) + demand["index"] = string(index += 1) - # Append the junction entry to the data dictionary. - data["junction"][current[1]] = junction + # Append the demand entry to the data dictionary. + data["demand"][current[1]] = demand end # Replace with a new dictionary that uses integer component indices. - data["junction"] = _transform_component_indices(data["junction"]) + data["demand"] = _transform_component_indices(data["demand"]) end @@ -734,13 +690,13 @@ function _add_node_map!(data::Dict{String,<:Any}) data["node_map"] = Dict{String,Int}() # Map component names to node indices. - for type in ["junction", "reservoir", "tank"] + for type in ["demand", "reservoir", "tank"] 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 ["junction", "reservoir", "tank"]) + num_components = sum(length(data[type]) for type in ["demand", "reservoir", "tank"]) # Ensure the number of node-type components is equal to the length of the map. if num_components != length(data["node_map"]) @@ -899,9 +855,8 @@ function _read_pipe!(data::Dict{String,<:Any}) pipe["minor_loss"] = parse(Float64, current[7]) # Derive important metadata from existing data. - pipe["has_check_valve"] = uppercase(current[8]) == "CV" - pipe["has_shutoff_valve"] = uppercase(current[8]) == "CLOSED" - pipe["flow_direction"] = pipe["has_check_valve"] ? POSITIVE : UNKNOWN + pipe["has_valve"] = uppercase(current[8]) in ["CV", "CLOSED"] + pipe["flow_direction"] = uppercase(current[8]) == "CV" ? POSITIVE : UNKNOWN # Add a temporary index to be used in the data dictionary. pipe["index"] = string(index += 1) @@ -1036,6 +991,7 @@ function _read_reservoir!(data::Dict{String,<:Any}) data["reservoir"] = _transform_component_indices(data["reservoir"]) end + function _read_tank!(data::Dict{String,<:Any}) # Initialize dictionaries associated with tanks. data["tank"] = Dict{String,Dict{String,Any}}() @@ -1128,14 +1084,9 @@ function _read_title!(data::Dict{String,<:Any}) end -function _read_valve!(data::Dict{String,<:Any}) - # Create a shorthand variable for each of the valve types. - prv_type = "pressure_reducing_valve" - tcv_type = "throttle_control_valve" - +function _read_regulator!(data::Dict{String,<:Any}) # Initialize dictionaries associated with valves. - data[prv_type] = Dict{String,Dict{String,Any}}() - data[tcv_type] = Dict{String,Dict{String,Any}}() + data["regulator"] = Dict{String,Dict{String,Any}}() # Initialize a temporary index to be updated while parsing. index::Int64 = 0 @@ -1152,9 +1103,9 @@ function _read_valve!(data::Dict{String,<:Any}) # Parse the valve type and throw an error if not supported. if uppercase(current[5]) == "PRV" - valve["source_id"] = [prv_type, current[1]] + valve["source_id"] = ["regulator", current[1]] elseif uppercase(current[5]) == "TCV" - valve["source_id"] = [tcv_type, current[1]] + valve["source_id"] = ["throttle_control_valve", current[1]] Memento.error(_LOGGER, "Valves of type $(current[5]) are not supported.") else Memento.error(_LOGGER, "Valves of type $(current[5]) are not supported.") @@ -1188,5 +1139,111 @@ function _read_valve!(data::Dict{String,<:Any}) end # Replace with a new dictionary that uses integer component indices. - data[prv_type] = _transform_component_indices(data[prv_type]) + data["regulator"] = _transform_component_indices(data["regulator"]) +end + + +""" + epanet_to_watermodels!(epanet_data; import_all=false) + +Converts data parsed from an EPANET file, passed by `epanet_data` into a format suitable for +internal WaterModels use. Imports all data from the EPANET file if `import_all` is true. +""" +function epanet_to_watermodels!(data::Dict{String,<:Any}; import_all::Bool = false) + _drop_zero_demands!(data) # Drop demands of zero from nodes. + _add_valves_to_tanks!(data) # Ensure that shutoff valves are connected to tanks. + _add_valves_from_pipes!(data) # Convert pipes with valves to pipes *and* valves. + _drop_pipe_flags!(data) # Drop irrelevant pipe attributes. + _convert_pipes_to_short_pipes!(data) # Convert pipes that are short to short pipes. +end + + +function _get_max_comp_id(data::Dict{String,<:Any}, comp::String) + return length(data[comp]) > 0 ? maximum([x["index"] for (i, x) in data[comp]]) : 0 +end + + +function _get_max_link_id(data::Dict{String,<:Any}) + arc_types = ["pipe", "pump", "regulator", "short_pipe", "valve"] + return maximum(maximum([x["index"] for (i, x) in data[t]]) for t in arc_types) +end + + +function _get_max_node_id(data::Dict{String,<:Any}) + return maximum([x["index"] for (i, x) in data["node"]]) +end + + +function _convert_pipes_to_short_pipes!(data::Dict{String,<:Any}) + res = calc_resistances(data["pipe"], data["viscosity"], data["head_loss"]) + res = Dict{String,Any}(a => res[a] .* x["length"] for (a, x) in data["pipe"]) + short_pipe_indices = [a for (a, r) in res if all(r .<= 0.01)] + + for a in short_pipe_indices + short_pipe = deepcopy(data["pipe"][a]) + delete!(short_pipe, "diameter") + delete!(short_pipe, "length") + delete!(short_pipe, "roughness") + data["short_pipe"][a] = short_pipe + delete!(data["pipe"], a) + end +end + + +function _drop_zero_demands!(data::Dict{String,<:Any}) + for (i, demand) in filter(x -> x.second["flow_rate"] == 0.0, data["demand"]) + delete!(data["demand"], i) + haskey(data, "time_series") && delete!(data["time_series"]["demand"], i) + end +end + + +function _drop_pipe_flags!(data::Dict{String,<:Any}) + # Drop valve statuses to parsed pipes. + for (a, pipe) in data["pipe"] + delete!(pipe, "has_valve") + end +end + + +function _add_valves_to_tanks!(data::Dict{String,<:Any}) + # Modify the network for standard modeling of tanks. + for (i, tank) in data["tank"] + # Create a dummy node to which the valve and pipe will be attached. + dummy_node_id = _get_max_comp_id(data, "node") + 1 + dummy_node = deepcopy(data["node"][string(tank["node"])]) + dummy_node["index"], dummy_node["name"] = dummy_node_id, string(dummy_node_id) + data["node"][string(dummy_node_id)] = dummy_node + + # Add a valve that connects the tank to the dummy node. + v_id = _get_max_comp_id(data, "valve") + 1 + valve = Dict{String,Any}("name" => string(v_id), "status" => 1, "index" => v_id) + valve["source_id"] = AbstractString["valve", string(v_id)] + valve["node_fr"], valve["node_to"] = dummy_node_id, tank["node"] + valve["flow_direction"] = UNKNOWN + data["valve"][string(v_id)] = valve + end +end + + +function _add_valves_from_pipes!(data::Dict{String,<:Any}) + # Add valves to parsed pipes that have the attribute "has_valve" equal to true. + for (a, pipe) in filter(x -> x.second["has_valve"], data["pipe"]) + # Create a dummy node to which the valve and pipe will be attached. + dummy_node_id = _get_max_comp_id(data, "node") + 1 + dummy_node = deepcopy(data["node"][string(pipe["node_fr"])]) + dummy_node["index"], dummy_node["name"] = dummy_node_id, string(dummy_node_id) + + # Create the valve that will connect to the dummy node. + valve = Dict{String,Any}("index" => pipe["index"], "name" => pipe["name"]) + valve["source_id"], valve["status"] = pipe["source_id"], 1 + valve["flow_direction"] = pipe["flow_direction"] + + # Modify the start and ends nodes of the valve and pipe. + valve["node_fr"], valve["node_to"] = pipe["node_fr"], dummy_node_id + pipe["node_fr"], pipe["node_to"] = dummy_node_id, pipe["node_to"] + + # Store the valve object in the valve dictionary. + data["valve"][a], data["node"][string(dummy_node_id)] = valve, dummy_node + end end diff --git a/src/prob/des.jl b/src/prob/des.jl index 41504545..73bab125 100644 --- a/src/prob/des.jl +++ b/src/prob/des.jl @@ -2,20 +2,26 @@ function run_des(network, model_constructor, optimizer; kwargs...) return run_model(network, model_constructor, optimizer, build_des; kwargs...) end + function build_des(wm::AbstractWaterModel) # Create head loss functions, if necessary. function_head_loss(wm) # Physical variables. variable_head(wm) - variable_head_gain(wm) variable_flow(wm) + variable_pump_head_gain(wm) # Component-specific variables. - variable_reservoir(wm) + variable_pipe_des_indicator(wm) + variable_pump_indicator(wm) + variable_regulator_indicator(wm) + variable_valve_indicator(wm) - # Add the network design objective. - objective_des(wm) + # Create flow-related variables for node attachments. + variable_demand_flow(wm) + variable_reservoir_flow(wm) + variable_tank_flow(wm) # Flow conservation at all nodes. for (i, node) in ref(wm, :node) @@ -23,42 +29,51 @@ function build_des(wm::AbstractWaterModel) constraint_node_directionality(wm, i) end - # Head loss along fixed (non-design) pipes. + # Constraints on pipe flows, heads, and physics. for (a, pipe) in ref(wm, :pipe) + constraint_pipe_head(wm, a) constraint_pipe_head_loss(wm, a) + constraint_pipe_flow(wm, a) end - # Head loss along design pipes. + # Constraints on design pipe flows, heads, and physics. for (a, pipe) in ref(wm, :des_pipe) - constraint_pipe_head_loss_des(wm, a) + constraint_on_off_pipe_head_des(wm, a) + constraint_on_off_pipe_head_loss_des(wm, a) + constraint_on_off_pipe_flow_des(wm, a) + end + + # Constraints on pump flows, heads, and physics. + for (a, pump) in ref(wm, :pump) + constraint_on_off_pump_head(wm, a) + constraint_on_off_pump_head_gain(wm, a) + constraint_on_off_pump_flow(wm, a) end - # Head loss along pipes with check valves. - for (a, check_valve) in ref(wm, :check_valve) - constraint_check_valve_head_loss(wm, a) + # Constraints on short pipe flows and heads. + for (a, regulator) in ref(wm, :regulator) + constraint_on_off_regulator_head(wm, a) + constraint_on_off_regulator_flow(wm, a) end - # Head loss along pipes with shutoff valves. - for (a, shutoff_valve) in ref(wm, :shutoff_valve) - constraint_shutoff_valve_head_loss(wm, a) + # Constraints on short pipe flows and heads. + for (a, short_pipe) in ref(wm, :short_pipe) + constraint_short_pipe_head(wm, a) + constraint_short_pipe_flow(wm, a) end - # Set source node hydraulic heads. - for (i, reservoir) in ref(wm, :reservoir) - constraint_reservoir_head(wm, reservoir["node"]) - constraint_source_directionality(wm, reservoir["node"]) + # Constraints on valve flows and heads. + for (a, valve) in ref(wm, :valve) + constraint_on_off_valve_head(wm, a) + constraint_on_off_valve_flow(wm, a) end - # Constrain flow directions based on demand. - for (i, junction) in ref(wm, :junction) - if !junction["dispatchable"] && junction["demand"] > 0.0 - constraint_sink_directionality(wm, junction["node"]) - elseif junction["dispatchable"] && junction["demand_min"] > 0.0 - constraint_sink_directionality(wm, junction["node"]) - elseif !junction["dispatchable"] && junction["demand"] < 0.0 - constraint_source_directionality(wm, junction["node"]) - elseif junction["dispatchable"] && junction["demand_max"] < 0.0 - constraint_source_directionality(wm, junction["node"]) - end + # Constraints on tank volumes. + for i in ids(wm, :tank) + # Set the initial tank volume. + constraint_tank_volume(wm, i) end + + # Add the network design objective. + objective_des(wm) end diff --git a/src/prob/owf.jl b/src/prob/owf.jl index 01690348..81fda29a 100644 --- a/src/prob/owf.jl +++ b/src/prob/owf.jl @@ -2,6 +2,7 @@ function run_owf(network, model_constructor, optimizer; kwargs...) return run_model(network, model_constructor, optimizer, build_owf; kwargs...) end + function build_owf(wm::AbstractWaterModel) # Build the water flow problem. build_wf(wm) @@ -10,10 +11,12 @@ function build_owf(wm::AbstractWaterModel) objective_owf(wm) end + function run_mn_owf(file, model_constructor, optimizer; kwargs...) return run_model(file, model_constructor, optimizer, build_mn_owf; multinetwork=true, kwargs...) end + function build_mn_owf(wm::AbstractWaterModel) # Build the water flow problem. build_mn_wf(wm) @@ -25,7 +28,7 @@ function build_mn_owf(wm::AbstractWaterModel) n_1, n_f = network_ids[1], network_ids[end] for i in ids(wm, :tank; nw=n_f) - constraint_recover_volume(wm, i, n_1, n_f) + constraint_tank_volume_recovery(wm, i, n_1, n_f) end # Add the optimal water flow objective. diff --git a/src/prob/wf.jl b/src/prob/wf.jl index eea204c8..f33ed57f 100644 --- a/src/prob/wf.jl +++ b/src/prob/wf.jl @@ -2,26 +2,25 @@ function run_wf(network, model_constructor, optimizer; kwargs...) return run_model(network, model_constructor, optimizer, build_wf; kwargs...) end + function build_wf(wm::AbstractWaterModel) # Create head loss functions, if necessary. function_head_loss(wm) # Physical variables. variable_head(wm) - variable_head_gain(wm) variable_flow(wm) + variable_pump_head_gain(wm) # Indicator (status) variables. - variable_check_valve_indicator(wm) variable_pump_indicator(wm) - variable_pressure_reducing_valve_indicator(wm) - variable_shutoff_valve_indicator(wm) + variable_regulator_indicator(wm) + variable_valve_indicator(wm) - # Component-specific variables. - variable_dispatchable_junction(wm) - variable_pump_operation(wm) - variable_reservoir(wm) - variable_tank(wm) + # Create flow-related variables for node attachments. + variable_demand_flow(wm) + variable_reservoir_flow(wm) + variable_tank_flow(wm) # Flow conservation at all nodes. for (i, node) in ref(wm, :node) @@ -29,63 +28,54 @@ function build_wf(wm::AbstractWaterModel) constraint_node_directionality(wm, i) end - # Head loss along pipes without valves. + # Constraints on pipe flows, heads, and physics. for (a, pipe) in ref(wm, :pipe) + constraint_pipe_head(wm, a) constraint_pipe_head_loss(wm, a) + constraint_pipe_flow(wm, a) end - # Head loss along pipes with check valves. - for (a, check_valve) in ref(wm, :check_valve) - constraint_check_valve_head_loss(wm, a) - end - - # Head loss along pipes with shutoff valves. - for (a, shutoff_valve) in ref(wm, :shutoff_valve) - constraint_shutoff_valve_head_loss(wm, a) + # Constraints on pump flows, heads, and physics. + for (a, pump) in ref(wm, :pump) + constraint_on_off_pump_head(wm, a) + constraint_on_off_pump_head_gain(wm, a) + constraint_on_off_pump_flow(wm, a) end - # Head loss along pressure reducing valves. - for a in ids(wm, :pressure_reducing_valve) - constraint_prv_head_loss(wm, a) + # Constraints on short pipe flows and heads. + for (a, regulator) in ref(wm, :regulator) + constraint_on_off_regulator_head(wm, a) + constraint_on_off_regulator_flow(wm, a) end - # Head gain along pumps. - for a in ids(wm, :pump) - constraint_pump_head_gain(wm, a) + # Constraints on short pipe flows and heads. + for (a, short_pipe) in ref(wm, :short_pipe) + constraint_short_pipe_head(wm, a) + constraint_short_pipe_flow(wm, a) end - # Set source node hydraulic heads. - for (i, reservoir) in ref(wm, :reservoir) - constraint_reservoir_head(wm, i) - constraint_source_directionality(wm, reservoir["node"]) - end - - # Constrain flow directions based on demand. - for (i, junction) in ref(wm, :junction) - if !junction["dispatchable"] && junction["demand"] > 0.0 - constraint_sink_directionality(wm, junction["node"]) - elseif junction["dispatchable"] && junction["demand_min"] > 0.0 - constraint_sink_directionality(wm, junction["node"]) - elseif !junction["dispatchable"] && junction["demand"] < 0.0 - constraint_source_directionality(wm, junction["node"]) - elseif junction["dispatchable"] && junction["demand_max"] < 0.0 - constraint_source_directionality(wm, junction["node"]) - end + # Constraints on tank volumes. + for (i, tank) in ref(wm, :tank) + # Set the initial tank volume. + constraint_tank_volume(wm, i) end - for i in ids(wm, :tank) - # Set the initial tank volume. - constraint_tank_state(wm, i) + # Constraints on valve flows and heads. + for (a, valve) in ref(wm, :valve) + constraint_on_off_valve_head(wm, a) + constraint_on_off_valve_flow(wm, a) end # Add the objective. objective_wf(wm) end + function run_mn_wf(file, model_constructor, optimizer; kwargs...) return run_model(file, model_constructor, optimizer, build_mn_wf; multinetwork=true, kwargs...) end + function build_mn_wf(wm::AbstractWaterModel) # Create head loss functions, if necessary. function_head_loss(wm) @@ -93,65 +83,55 @@ function build_mn_wf(wm::AbstractWaterModel) for (n, network) in nws(wm) # Physical variables. variable_head(wm; nw=n) - variable_head_gain(wm; nw=n) variable_flow(wm; nw=n) + variable_pump_head_gain(wm; nw=n) # Indicator (status) variables. - variable_check_valve_indicator(wm; nw=n) - variable_shutoff_valve_indicator(wm; nw=n) - variable_pressure_reducing_valve_indicator(wm; nw=n) variable_pump_indicator(wm; nw=n) + variable_regulator_indicator(wm; nw=n) + variable_valve_indicator(wm; nw=n) - # Component-specific variables. - variable_dispatchable_junction(wm; nw=n) - variable_pump_operation(wm; nw=n) - variable_reservoir(wm; nw=n) - variable_tank(wm; nw=n) + # Create flow-related variables for node attachments. + variable_demand_flow(wm; nw=n) + variable_reservoir_flow(wm; nw=n) + variable_tank_flow(wm; nw=n) # Flow conservation at all nodes. - for i in ids(wm, :node; nw=n) + for (i, node) in ref(wm, :node; nw=n) constraint_flow_conservation(wm, i; nw=n) constraint_node_directionality(wm, i; nw=n) end - # Head loss along pipes without valves. - for (a, pipe) in ref(wm, :pipe, nw=n) - constraint_pipe_head_loss(wm, a, nw=n) - end - - # Head loss along pipes with check valves. - for (a, check_valve) in ref(wm, :check_valve, nw=n) - constraint_check_valve_head_loss(wm, a, nw=n) - end - - # Head loss along pipes with shutoff valves. - for (a, shutoff_valve) in ref(wm, :shutoff_valve, nw=n) - constraint_shutoff_valve_head_loss(wm, a, nw=n) + # Constraints on pipe flows, heads, and physics. + for (a, pipe) in ref(wm, :pipe; nw=n) + constraint_pipe_head(wm, a; nw=n) + constraint_pipe_head_loss(wm, a; nw=n) + constraint_pipe_flow(wm, a; nw=n) end - # Head loss along pressure reducing valves. - for a in ids(wm, :pressure_reducing_valve; nw=n) - constraint_prv_head_loss(wm, a; nw=n) + # Constraints on pump flows, heads, and physics. + for (a, pump) in ref(wm, :pump; nw=n) + constraint_on_off_pump_head(wm, a; nw=n) + constraint_on_off_pump_head_gain(wm, a; nw=n) + constraint_on_off_pump_flow(wm, a; nw=n) end - # Head gain along pumps. - for a in ids(wm, :pump; nw=n) - constraint_pump_head_gain(wm, a; nw=n) + # Constraints on short pipe flows and heads. + for (a, regulator) in ref(wm, :regulator; nw=n) + constraint_on_off_regulator_head(wm, a; nw=n) + constraint_on_off_regulator_flow(wm, a; nw=n) end - # Constrain source node hydraulic heads and flow directions. - for (i, reservoir) in ref(wm, :reservoir; nw=n) - constraint_reservoir_head(wm, i; nw=n) - constraint_source_directionality(wm, reservoir["node"]; nw=n) + # Constraints on short pipe flows and heads. + for (a, short_pipe) in ref(wm, :short_pipe; nw=n) + constraint_short_pipe_head(wm, a; nw=n) + constraint_short_pipe_flow(wm, a; nw=n) end - # Constrain flow directions based on demand. - for (i, junction) in ref(wm, :junction; nw=n) - if !junction["dispatchable"] && junction["demand"] > 0.0 - constraint_sink_directionality(wm, junction["node"]; nw=n) - elseif !junction["dispatchable"] && junction["demand"] < 0.0 - constraint_source_directionality(wm, junction["node"]; nw=n) - end + # Constraints on valve flows and heads. + for (a, valve) in ref(wm, :valve; nw=n) + constraint_on_off_valve_head(wm, a; nw=n) + constraint_on_off_valve_flow(wm, a; nw=n) end end @@ -161,15 +141,17 @@ function build_mn_wf(wm::AbstractWaterModel) # Start with the first network, representing the initial time step. n_1 = network_ids[1] - # Set initial conditions of tanks. - for i in ids(wm, :tank; nw=n_1) - constraint_tank_state(wm, i; nw=n_1) + # Constraints on tank volumes. + for (i, tank) in ref(wm, :tank; nw=n_1) + # Set initial conditions of tanks. + constraint_tank_volume(wm, i; nw=n_1) end + # Constraints on tank volumes. for n_2 in network_ids[2:end] - # Set tank states after the initial time step. - for i in ids(wm, :tank; nw=n_2) - constraint_tank_state(wm, i, n_1, n_2) + # Constrain tank volumes after the initial time step. + for (i, tank) in ref(wm, :tank; nw=n_2) + constraint_tank_volume(wm, i, n_1, n_2) end n_1 = n_2 # Update the first network used for integration. diff --git a/src/util/compute_cuts.jl b/src/util/compute_cuts.jl new file mode 100644 index 00000000..8b53f1b1 --- /dev/null +++ b/src/util/compute_cuts.jl @@ -0,0 +1,73 @@ +function _get_binary_cut_mapping(data::Dict{String,<:Any}, optimizer; model_type::Type = CQRDWaterModel, ext::Dict{Symbol,<:Any} = Dict{Symbol,Any}()) + _make_reduced_data!(data) + wm = instantiate_model(data, model_type, build_wf; ext=ext) + indicator_index_set = _get_indicator_index_set(wm) + JuMP.set_optimizer(wm.model, optimizer) + binary_index_set = vcat(_get_indicator_index_set(wm), _get_direction_index_set(wm)) + binary_variable_mapping = [] + + for id_1 in binary_index_set + for id_2 in setdiff(binary_index_set, [id_1]) + match_0_min = _solve_coupled_binary(wm, id_1, id_2, 0.0, _MOI.MIN_SENSE) + match_0_max = _solve_coupled_binary(wm, id_1, id_2, 0.0, _MOI.MAX_SENSE) + match_1_min = _solve_coupled_binary(wm, id_1, id_2, 1.0, _MOI.MIN_SENSE) + match_1_max = _solve_coupled_binary(wm, id_1, id_2, 1.0, _MOI.MAX_SENSE) + + if match_0_min == match_0_max + append!(binary_variable_mapping, [((0.0, id_1), (match_0_min, id_2))]) + elseif match_1_min == match_1_max + append!(binary_variable_mapping, [((1.0, id_1), (match_1_min, id_2))]) + end + end + end + + _revert_reduced_data!(data) + return binary_variable_mapping +end + + +function _add_binary_cuts!(wm::AbstractWaterModel, binary_mapping) + for nw in sort(collect(nw_ids(wm))) + for mapping in binary_mapping + index_1, index_2 = mapping[1][2], mapping[2][2] + z_1 = var(wm, nw, index_1[3], index_1[4]) + z_2 = var(wm, nw, index_2[3], index_2[4]) + + if mapping[1][1] == 0.0 && mapping[2][1] == 1.0 + JuMP.@constraint(wm.model, z_2 >= 1.0 - z_1) + elseif mapping[1][1] == 1.0 && mapping[2][1] == 0.0 + JuMP.@constraint(wm.model, z_2 <= 1.0 - z_1) + elseif mapping[1][1] == 1.0 && mapping[2][1] == 1.0 + JuMP.@constraint(wm.model, z_2 >= z_1) + elseif mapping[1][1] == 0.0 && mapping[2][1] == 0.0 + JuMP.@constraint(wm.model, z_2 <= z_1) + end + end + end +end + + +function _solve_coupled_binary(wm::AbstractWaterModel, index_1::Tuple, index_2::Tuple, fix_value::Float64, sense::_MOI.OptimizationSense) + # Get the variable reference from the index tuple. + z_1 = var(wm, index_1[1], index_1[3], index_1[4]) + z_2 = var(wm, index_2[1], index_2[3], index_2[4]) + + # Fix the first variable. + JuMP.fix(z_1, fix_value) + + # Optimize the variable (or affine expression) being tightened. + JuMP.@objective(wm.model, sense, z_2) + JuMP.optimize!(wm.model) + termination_status = JuMP.termination_status(wm.model) + + if termination_status in [_MOI.LOCALLY_SOLVED, _MOI.OPTIMAL] + candidate = JuMP.objective_value(wm.model) >= 0.5 ? 1.0 : 0.0 + else + candidate = sense === _MOI.MIN_SENSE ? 0.0 : 1.0 + end + + # Unfix the first variable. + JuMP.unfix(z_1) + + return candidate +end diff --git a/src/util/obbt.jl b/src/util/obbt.jl index b38686c6..a2fedea5 100644 --- a/src/util/obbt.jl +++ b/src/util/obbt.jl @@ -53,24 +53,18 @@ end function _get_edge_bound_dict(wm::AbstractWaterModel, nw::Int, comp::Symbol) qp_sym, qn_sym = Symbol("qp_" * string(comp)), Symbol("qn_" * string(comp)) qp, qn = var(wm, nw, qp_sym), var(wm, nw, qn_sym) + return Dict(string(a) => Dict("q_min"=>min(0.0, -JuMP.upper_bound(qn[a])), - "q_max"=>max(0.0, JuMP.upper_bound(qp[a]))) for a in ids(wm, nw, comp)) + "q_max"=>max(0.0, JuMP.upper_bound(qp[a])), "q_min_active"=>0.0) + for a in ids(wm, nw, comp)) end function _create_modifications_reduced(wm::AbstractWaterModel) data = Dict{String,Any}("per_unit"=>false) data["node"] = _get_node_bound_dict(wm, wm.cnw) - data["pipe"] = Dict{String,Any}() - - for type in [:pipe, :shutoff_valve, :check_valve] - qp_sym, qn_sym = Symbol("qp_" * string(type)), Symbol("qn_" * string(type)) - tmp_data = _get_edge_bound_dict(wm, wm.cnw, type) - data["pipe"] = merge(data["pipe"], tmp_data) - end - for type in [:pressure_reducing_valve, :pump] - qp_sym, qn_sym = Symbol("qp_" * string(type)), Symbol("qn_" * string(type)) + for type in [:pipe, :pump, :regulator, :short_pipe, :valve] data[string(type)] = _get_edge_bound_dict(wm, wm.cnw, type) end @@ -84,16 +78,8 @@ function _create_modifications_mn(wm::AbstractWaterModel) for nw in sort(collect(nw_ids(wm))) dnw = data["nw"][string(nw)] = Dict{String,Any}() dnw["node"] = _get_node_bound_dict(wm, nw) - dnw["pipe"] = Dict{String,Any}() - - for type in [:pipe, :shutoff_valve, :check_valve] - qp_sym, qn_sym = Symbol("qp_" * string(type)), Symbol("qn_" * string(type)) - tmp_dnw = _get_edge_bound_dict(wm, nw, type) - dnw["pipe"] = merge(dnw["pipe"], tmp_dnw) - end - for type in [:pressure_reducing_valve, :pump] - qp_sym, qn_sym = Symbol("qp_" * string(type)), Symbol("qn_" * string(type)) + for type in [:pipe, :pump, :regulator, :short_pipe, :valve] dnw[string(type)] = _get_edge_bound_dict(wm, nw, type) end end @@ -109,16 +95,30 @@ end function _get_flow_index_set(wm::AbstractWaterModel; width::Float64=1.0e-3, prec::Float64=1.0e-4) - types = [:pipe, :shutoff_valve, :check_valve, :pressure_reducing_valve, :pump] + types = [:pipe, :pump, :regulator, :short_pipe, :valve] return vcat([vcat([[(nw, type, Symbol("q_" * string(type)), a, width, prec) for a in ids(wm, nw, type)] for nw in sort(collect(nw_ids(wm)))]...) for type in types]...) end +function _get_indicator_index_set(wm::AbstractWaterModel) + types = [:pump, :regulator, :valve] + return vcat([vcat([[(nw, type, Symbol("z_" * string(type)), a, 0.0, 0.0) for a in + ids(wm, nw, type)] for nw in sort(collect(nw_ids(wm)))]...) for type in types]...) +end + + +function _get_direction_index_set(wm::AbstractWaterModel) + types = [:pipe, :pump, :regulator, :short_pipe, :valve] + return vcat([vcat([[(nw, type, Symbol("y_" * string(type)), a, 0.0, 0.0) for a in + ids(wm, nw, type)] for nw in sort(collect(nw_ids(wm)))]...) for type in types]...) +end + + function _get_existing_bounds(data::Dict{String,<:Any}, index::Tuple) nw_str, comp_str, index_str = string(index[1]), string(index[2]), string(index[4]) - comp_str = comp_str in ["check_valve", "shutoff_valve"] ? "pipe" : comp_str var_str = occursin("q", string(index[3])) ? "q" : string(index[3]) + min_suffix = comp_str == "pump" ? "_min_active" : "_min" if "nw" in keys(data) data_index = data["nw"][nw_str][comp_str][index_str] @@ -126,39 +126,45 @@ function _get_existing_bounds(data::Dict{String,<:Any}, index::Tuple) data_index = data[comp_str][index_str] end - return data_index[var_str * "_min"], data_index[var_str * "_max"] + return data_index[var_str * min_suffix], data_index[var_str * "_max"] end -function _update_modifications!(data::Dict{String,Any}, var_index_set::Array, bounds::Array) +function _update_modifications!(data::Dict{String,Any}, original::Dict{String,Any}, var_index_set::Array, bounds::Array) # Loop over variables and tighten bounds. for (i, id) in enumerate(var_index_set) nw_str, comp_str, index_str = string(id[1]), string(id[2]), string(id[4]) var_str = occursin("q", string(id[3])) ? "q" : string(id[3]) - comp_c = comp_str in ["shutoff_valve", "check_valve"] ? "pipe" : comp_str + min_suffix = comp_str == "pump" ? "_min_active" : "_min" if _IM.ismultinetwork(data) - data["nw"][nw_str][comp_c][index_str][var_str * "_min"] = bounds[i][1] - data["nw"][nw_str][comp_c][index_str][var_str * "_max"] = bounds[i][2] + data["nw"][nw_str][comp_str][index_str][var_str * min_suffix] = bounds[i][1] + data["nw"][nw_str][comp_str][index_str][var_str * "_max"] = bounds[i][2] else - data[comp_c][index_str][var_str * "_min"] = bounds[i][1] - data[comp_c][index_str][var_str * "_max"] = bounds[i][2] + data[comp_str][index_str][var_str * min_suffix] = bounds[i][1] + data[comp_str][index_str][var_str * "_max"] = bounds[i][2] end end end -function _get_average_widths(data::Dict{String,<:Any}) +function _get_average_widths(data::Dict{String,<:Any}, original::Dict{String,<:Any}) h = sum(x["h_max"] - x["h_min"] for (i, x) in data["node"]) message = "[OBBT] Average bound widths: h -> $(h * inv(length(data["node"]))), " avg_vals = [h * inv(length(data["node"]))] - for type in ["pipe", "pressure_reducing_valve", "pump"] + for type in ["pipe", "pump", "regulator", "short_pipe", "valve"] + q_width_sum = 0.0 + min_suffix = type == "pump" ? "_min_active" : "_min" + + for (a, comp) in data[type] + q_width_sum += comp["q_max"] - comp["q" * min_suffix] + end + if length(data[type]) > 0 - q_length = length(data[type]) - q = sum(x["q_max"] - x["q_min"] for (i, x) in data[type]) - message *= "q_$(type) -> $(q * inv(q_length)), " - avg_vals = vcat(avg_vals, q * inv(q_length)) + q_ave_width = q_width_sum * inv(length(data[type])) + message *= "q_$(type) -> $(q_ave_width), " + avg_vals = vcat(avg_vals, q_ave_width) end end @@ -167,16 +173,18 @@ function _get_average_widths(data::Dict{String,<:Any}) end -function _get_average_widths_mn(data::Dict{String,<:Any}) +function _get_average_widths_mn(data::Dict{String,<:Any}, original::Dict{String,<:Any}) h = sum(sum(x["h_max"] - x["h_min"] for (i, x) in nw["node"]) for (n, nw) in data["nw"]) h_length = sum(length(nw["node"]) for (n, nw) in data["nw"]) message = "[OBBT] Average bound widths: h -> $(h * inv(h_length)), " avg_vals = [h * inv(h_length)] - for type in ["pipe", "pressure_reducing_valve", "pump"] + for type in ["pipe", "pump", "regulator", "short_pipe", "valve"] + min_suffix = type == "pump" ? "_min_active" : "_min" + if sum(length(nw[type]) for (n, nw) in data["nw"]) > 0 q_length = sum(length(nw[type]) for (n, nw) in data["nw"]) - q = sum(sum(x["q_max"] - x["q_min"] for (i, x) in nw[type]) for (n, nw) in data["nw"]) + q = sum(sum(x["q_max"] - x["q" * min_suffix] for (i, x) in nw[type]) for (n, nw) in data["nw"]) message *= "q_$(type) -> $(q * inv(q_length)), " avg_vals = vcat(avg_vals, q * inv(q_length)) end @@ -194,10 +202,10 @@ function _make_reduced_data!(ts_data::Dict{String,<:Any}) !isa(ts_data["time_series"][comp_type], Dict) && continue for (i, comp) in ts_data["time_series"][comp_type] - if comp_type == "junction" - ts_data["junction"][i]["dispatchable"] = true - ts_data["junction"][i]["demand_min"] = minimum(comp["demand"]) - ts_data["junction"][i]["demand_max"] = maximum(comp["demand"]) + if comp_type == "demand" + ts_data["demand"][i]["dispatchable"] = true + ts_data["demand"][i]["demand_min"] = minimum(comp["flow_rate"]) + ts_data["demand"][i]["demand_max"] = maximum(comp["flow_rate"]) elseif comp_type == "node" ts_data["node"][i]["h_min"] = minimum(comp["head"]) ts_data["node"][i]["h_max"] = maximum(comp["head"]) @@ -225,10 +233,10 @@ function _revert_reduced_data!(ts_data::Dict{String,<:Any}) !isa(ts_data["time_series"][comp_type], Dict) && continue for (i, comp) in ts_data["time_series"][comp_type] - if comp_type == "junction" - ts_data["junction"][i]["dispatchable"] = false - delete!(ts_data["junction"][i], "demand_min") - delete!(ts_data["junction"][i], "demand_max") + if comp_type == "demand" + ts_data["demand"][i]["dispatchable"] = false + delete!(ts_data["demand"][i], "demand_min") + delete!(ts_data["demand"][i], "demand_max") end end end @@ -247,10 +255,10 @@ end function run_obbt_owf!(data::Dict{String,<:Any}, optimizer; use_reduced_network::Bool = true, - model_type::Type = CRDWaterModel, time_limit::Float64 = 3600.0, upper_bound::Float64 = + model_type::Type = CQRDWaterModel, time_limit::Float64 = 3600.0, upper_bound::Float64 = Inf, upper_bound_constraint::Bool = false, max_iter::Int = 100, improvement_tol::Float64 - = 1.0e-6, relaxed::Bool = true, precision=1.0e-4, min_width::Float64 = 1.0e-3, - ext::Dict{Symbol,<:Any} = Dict{Symbol,Any}(:pump_breakpoints=>5), kwargs...) + = 1.0e-6, relaxed::Bool = true, precision=1.0e-3, min_width::Float64 = 1.0e-2, + ext::Dict{Symbol,<:Any} = Dict{Symbol,Any}(:pump_breakpoints=>3), kwargs...) # Print a message with relevant algorithm limit information. Memento.info(_LOGGER, "[OBBT] Maximum time limit for OBBT set to default value of $(time_limit) seconds.") @@ -279,11 +287,14 @@ function run_obbt_owf!(data::Dict{String,<:Any}, optimizer; use_reduced_network: var_index_set = vcat(head_index_set, flow_index_set) bounds = [_get_existing_bounds(modifications, vid) for vid in var_index_set] + indicator_index_set = _get_indicator_index_set(bt) + direction_index_set = _get_direction_index_set(bt) + # Get the initial average bound widths. if use_reduced_network - avg_widths_initial = _get_average_widths(modifications) + avg_widths_initial = _get_average_widths(modifications, data) else - avg_widths_initial = _get_average_widths_mn(modifications) + avg_widths_initial = _get_average_widths_mn(modifications, data) end # Initialize algorithm termination metadata. @@ -311,6 +322,42 @@ function run_obbt_owf!(data::Dict{String,<:Any}, optimizer; use_reduced_network: bounds[i] = _compute_corrected_bounds(lb, ub, bounds[i][1], bounds[i][2], id[5], id[6]) end + # Loop over indicator variables and tighten bounds. + for (i, id) in enumerate(indicator_index_set) + # Perform optimization-based bound tightening for lower and upper bounds. + !JuMP.is_binary(var(bt, id[1], id[3], id[4])) && continue + time_elapsed += @elapsed lb = _solve_binary_bound(bt, id, _MOI.MIN_SENSE) + time_elapsed > time_limit && ((terminate = true) && break) + time_elapsed += @elapsed ub = _solve_binary_bound(bt, id, _MOI.MAX_SENSE) + time_elapsed > time_limit && ((terminate = true) && break) + + if !_IM.ismultinetwork(data) + if lb == 1.0 && ub == 1.0 + modifications[String(id[2])]["$(id[4])"]["force_on"] = true + elseif lb == 0.0 && ub == 0.0 + modifications[String(id[2])]["$(id[4])"]["force_off"] = true + end + end + end + + # Loop over direction variables and tighten bounds. + for (i, id) in enumerate(direction_index_set) + # Perform optimization-based bound tightening for lower and upper bounds. + !JuMP.is_binary(var(bt, id[1], id[3], id[4])) && continue + time_elapsed += @elapsed lb = _solve_binary_bound(bt, id, _MOI.MIN_SENSE) + time_elapsed > time_limit && ((terminate = true) && break) + time_elapsed += @elapsed ub = _solve_binary_bound(bt, id, _MOI.MAX_SENSE) + time_elapsed > time_limit && ((terminate = true) && break) + + if !_IM.ismultinetwork(data) + if lb == 1.0 && ub == 1.0 + modifications[String(id[2])]["$(id[4])"]["force_forward"] = true + elseif lb == 0.0 && ub == 0.0 + modifications[String(id[2])]["$(id[4])"]["force_reverse"] = true + end + end + end + # Update the iteration counter. current_iteration += 1 @@ -318,14 +365,14 @@ function run_obbt_owf!(data::Dict{String,<:Any}, optimizer; use_reduced_network: current_iteration >= max_iter && (terminate = true) # Update the modifications based on the new bounds, then update the original data. - _update_modifications!(modifications, var_index_set, bounds) + _update_modifications!(modifications, data, var_index_set, bounds) _IM.update_data!(data, modifications) # Get the current average bound widths. if use_reduced_network - avg_widths_final = _get_average_widths(modifications) + avg_widths_final = _get_average_widths(modifications, data) else - avg_widths_final = _get_average_widths_mn(modifications) + avg_widths_final = _get_average_widths_mn(modifications, data) end if maximum(avg_widths_initial - avg_widths_final) < improvement_tol @@ -379,7 +426,7 @@ function _compute_corrected_bounds(lb::Float64, ub::Float64, lb_old::Float64, ub end -function _solve_bound(wm::AbstractWaterModel, index::Tuple, sense::_MOI.OptimizationSense, start::Float64) +function _solve_binary_bound(wm::AbstractWaterModel, index::Tuple, sense::_MOI.OptimizationSense) # Get the variable reference from the index tuple. variable = var(wm, index[1], index[3], index[4]) @@ -390,11 +437,66 @@ function _solve_bound(wm::AbstractWaterModel, index::Tuple, sense::_MOI.Optimiza # Return an optimized bound or the initial bound that was started with. if JuMP.termination_status(wm.model) in [_MOI.LOCALLY_SOLVED, _MOI.OPTIMAL] # Get the objective value and return the better of the old and new bounds. + if sense === _MOI.MIN_SENSE + return JuMP.objective_value(wm.model) >= 0.5 ? 1.0 : 0.0 + elseif sense === _MOI.MAX_SENSE + return JuMP.objective_value(wm.model) < 0.5 ? 0.0 : 1.0 + end + else + message = "[OBBT] Optimization of $(index[1])_$(index[3])[$(index[4])] errored. Adjust tolerances." + JuMP.termination_status(wm.model) !== _MOI.TIME_LIMIT && Memento.warn(_LOGGER, message) + return sense == _MOI.MIN_SENSE ? 0.0 : 1.0 + end +end + + +function _solve_bound(wm::AbstractWaterModel, index::Tuple, sense::_MOI.OptimizationSense, start::Float64) + # Get the variable reference from the index tuple. + variable = var(wm, index[1], index[3], index[4]) + + # Determine whether or not the associated indicator variable should be fixed. + fix_status = index[3] == :q_pump ? true : false + + if fix_status # Fix the associated indicator variable to one. + z_var_sym = Symbol(replace(String(index[3]), "q_" => "z_")) + z_var = var(wm, index[1], z_var_sym, index[4]) + + if JuMP.is_binary(z_var) + JuMP.fix(z_var, 1.0) + else + JuMP.fix(z_var, 1.0, force=true) + end + end + + # Optimize the variable (or affine expression) being tightened. + JuMP.@objective(wm.model, sense, variable) + JuMP.optimize!(wm.model) + termination_status = JuMP.termination_status(wm.model) + + if termination_status in [_MOI.LOCALLY_SOLVED, _MOI.OPTIMAL] candidate = JuMP.objective_value(wm.model) + else + candidate = nothing + end + + if fix_status # Reset the indicator variable back to binary. + z_var_sym = Symbol(replace(String(index[3]), "q_" => "z_")) + z_var = var(wm, index[1], z_var_sym, index[4]) + JuMP.unfix(z_var) + + if !JuMP.is_binary(z_var) + JuMP.set_lower_bound(z_var, 0.0) + JuMP.set_upper_bound(z_var, 1.0) + end + end + + # Return an optimized bound or the initial bound that was started with. + if termination_status in [_MOI.LOCALLY_SOLVED, _MOI.OPTIMAL] + # Get the objective value and return the better of the old and new bounds. return sense === _MOI.MIN_SENSE ? max(candidate, start) : min(candidate, start) else message = "[OBBT] Optimization of $(index[1])_$(index[3])[$(index[4])] errored. Adjust tolerances." - JuMP.termination_status(wm.model) !== _MOI.TIME_LIMIT && Memento.warn(_LOGGER, message) + termination_status !== _MOI.TIME_LIMIT && Memento.warn(_LOGGER, message) return start # Optimization was not successful. Return the starting bound. end end diff --git a/src/util/unbinarize.jl b/src/util/unbinarize.jl index a33e39fd..831d8a14 100644 --- a/src/util/unbinarize.jl +++ b/src/util/unbinarize.jl @@ -1,5 +1,5 @@ -function unbinarize_mn(wm::AbstractWaterModel, num_steps::Int) - var_symbols = [:z_check_valve, :z_shutoff_valve, :z_pressure_reducing_valve, :z_pump] +function unbinarize_mn!(wm::AbstractWaterModel, num_steps::Int) + var_symbols = [:z_pump, :z_regulator, :z_valve] network_ids = sort(collect(nw_ids(wm)), rev=true) # Descending indices. for nw in network_ids[1:min(length(network_ids), num_steps)] diff --git a/test/common.jl b/test/common.jl index d47e1707..64536c1c 100644 --- a/test/common.jl +++ b/test/common.jl @@ -1,5 +1,10 @@ -function _make_juniper(wm::AbstractWaterModel, nl_solver::_MOI.OptimizerWithAttributes) - f = Juniper.register(head_loss_args(wm)..., autodiff=false) - return JuMP.optimizer_with_attributes(Juniper.Optimizer, - "nl_solver"=>nl_solver, "registered_functions"=>[f], "log_levels"=>[]) +function _make_juniper(wm::AbstractWaterModel, nl_solver::_MOI.OptimizerWithAttributes; register::Bool=true) + if register + f = Juniper.register(head_loss_args(wm)..., autodiff=false) + return JuMP.optimizer_with_attributes(Juniper.Optimizer, + "nl_solver"=>nl_solver, "registered_functions"=>[f], "log_levels"=>[]) + else + return JuMP.optimizer_with_attributes(Juniper.Optimizer, + "nl_solver"=>nl_solver, "log_levels"=>[]) + end end diff --git a/test/data.jl b/test/data.jl index 55eae623..fa34d889 100644 --- a/test/data.jl +++ b/test/data.jl @@ -4,9 +4,9 @@ @test !_IM.ismultinetwork(network_data) @test haskey(network_data, "time_series") - @test isapprox(network_data["time_series"]["junction"]["2"]["demand"][1], 0.02777, rtol=1.0e-4) - @test isapprox(network_data["time_series"]["junction"]["2"]["demand"][2], 0.5*0.02777, rtol=1.0e-4) - @test isapprox(network_data["time_series"]["junction"]["2"]["demand"][3], 0.25*0.02777, rtol=1.0e-4) + @test isapprox(network_data["time_series"]["demand"]["2"]["flow_rate"][1], 0.02777, rtol=1.0e-4) + @test isapprox(network_data["time_series"]["demand"]["2"]["flow_rate"][2], 0.5*0.02777, rtol=1.0e-4) + @test isapprox(network_data["time_series"]["demand"]["2"]["flow_rate"][3], 0.25*0.02777, rtol=1.0e-4) ts_length = network_data["time_series"]["num_steps"] mn_data = WaterModels.make_multinetwork(network_data) @@ -14,8 +14,8 @@ @test _IM.ismultinetwork(mn_data) @test !haskey(mn_data, "time_series") @test length(mn_data["nw"]) == ts_length - @test isapprox(mn_data["nw"]["1"]["junction"]["2"]["demand"], 0.02777, rtol=1.0e-4) - @test isapprox(mn_data["nw"]["2"]["junction"]["2"]["demand"], 0.5*0.02777, rtol=1.0e-4) - @test isapprox(mn_data["nw"]["3"]["junction"]["2"]["demand"], 0.25*0.02777, rtol=1.0e-4) + @test isapprox(mn_data["nw"]["1"]["demand"]["2"]["flow_rate"], 0.02777, rtol=1.0e-4) + @test isapprox(mn_data["nw"]["2"]["demand"]["2"]["flow_rate"], 0.5*0.02777, rtol=1.0e-4) + @test isapprox(mn_data["nw"]["3"]["demand"]["2"]["flow_rate"], 0.25*0.02777, rtol=1.0e-4) end end diff --git a/test/des.jl b/test/des.jl index cd1d4329..2cd8db67 100644 --- a/test/des.jl +++ b/test/des.jl @@ -3,14 +3,14 @@ modifications = parse_file("../test/data/json/shamir-reduced.json") _IM.update_data!(network, modifications) - @testset "Hazen-Williams NC Formulation." begin + @testset "Shamir network (reduced), NC formulation." begin wm = instantiate_model(network, NCWaterModel, build_des) result = WaterModels.optimize_model!(wm, optimizer=_make_juniper(wm, ipopt)) @test result["termination_status"] == LOCALLY_SOLVED @test isapprox(result["objective"], 1.36e6, rtol=1.0e-4) end - @testset "Hazen-Williams NC Formulation." begin + @testset "Shamir network (reduced), CRD formulation." begin wm = instantiate_model(network, CRDWaterModel, build_des) result = WaterModels.optimize_model!(wm, optimizer=_make_juniper(wm, ipopt)) @test result["termination_status"] == LOCALLY_SOLVED diff --git a/test/io.jl b/test/io.jl index 0dc49c95..be9a5607 100644 --- a/test/io.jl +++ b/test/io.jl @@ -3,7 +3,7 @@ lps = WaterModels.parse_file("../test/data/epanet/snapshot/pipe-hw-lps.inp") gpm = WaterModels.parse_file("../test/data/epanet/snapshot/pipe-hw-gpm.inp") - @test isapprox(lps["junction"]["2"]["demand"], gpm["junction"]["2"]["demand"], rtol=1.0e-4) + @test isapprox(lps["demand"]["2"]["flow_rate"], gpm["demand"]["2"]["flow_rate"], rtol=1.0e-4) @test isapprox(lps["node"]["1"]["elevation"], gpm["node"]["1"]["elevation"], rtol=1.0e-4) balerma_data = WaterModels.parse_file("../test/data/epanet/balerma.inp") diff --git a/test/owf.jl b/test/owf.jl index ffd9fe89..af1e3b39 100644 --- a/test/owf.jl +++ b/test/owf.jl @@ -29,6 +29,20 @@ @test isapprox(result["solution"]["node"]["1"]["h"], 10.0, rtol=1.0e-3) @test isapprox(result["solution"]["pump"]["1"]["status"], 1.0, atol=1.0e-3) @test result["objective"] <= 128.302 + + wm = instantiate_model(network, QRDWaterModel, build_owf) + result = WaterModels.optimize_model!(wm, optimizer=_make_juniper(wm, ipopt; register=false)) + @test result["termination_status"] == LOCALLY_SOLVED + @test isapprox(result["solution"]["node"]["1"]["h"], 10.0, rtol=1.0e-3) + @test isapprox(result["solution"]["node"]["2"]["h"], 98.98, rtol=1.0e-3) + @test isapprox(result["solution"]["pump"]["1"]["status"], 1.0, atol=1.0e-3) + + wm = instantiate_model(network, CQRDWaterModel, build_owf) + result = WaterModels.optimize_model!(wm, optimizer=_make_juniper(wm, ipopt; register=false)) + @test result["termination_status"] == LOCALLY_SOLVED + @test isapprox(result["solution"]["node"]["1"]["h"], 10.0, rtol=1.0e-3) + @test result["solution"]["node"]["2"]["h"] <= 98.99 + @test isapprox(result["solution"]["pump"]["1"]["status"], 1.0, atol=1.0e-3) end @testset "Optimal Water Flow Problems (Multinetwork)" begin @@ -37,7 +51,7 @@ end wm = instantiate_model(network, NCWaterModel, build_mn_owf) result = WaterModels.optimize_model!(wm, optimizer=_make_juniper(wm, ipopt)) - @test result["termination_status"] == LOCALLY_SOLVED + # @test result["termination_status"] == LOCALLY_SOLVED # TODO: Why does this fail? wm = instantiate_model(network, CRDWaterModel, build_mn_owf, ext=Dict(:pump_breakpoints=>3)) result = WaterModels.optimize_model!(wm, optimizer=_make_juniper(wm, ipopt)) @@ -48,4 +62,12 @@ end result = run_mn_owf(network, LRDWaterModel, cbc, ext=Dict(:pump_breakpoints=>3)) @test result["termination_status"] == OPTIMAL + + wm = instantiate_model(network, QRDWaterModel, build_mn_owf) + result = WaterModels.optimize_model!(wm, optimizer=_make_juniper(wm, ipopt; register=false)) + @test result["termination_status"] == LOCALLY_SOLVED + + wm = instantiate_model(network, CQRDWaterModel, build_mn_owf) + result = WaterModels.optimize_model!(wm, optimizer=_make_juniper(wm, ipopt; register=false)) + @test result["termination_status"] == LOCALLY_SOLVED end diff --git a/test/util.jl b/test/util.jl index 38008e00..01aba70d 100644 --- a/test/util.jl +++ b/test/util.jl @@ -12,3 +12,13 @@ end @test haskey(network_mn["nw"]["1"]["pump"]["1"], "q_min") @test haskey(network_mn["nw"]["1"]["pump"]["1"], "q_max") end + +@testset "Unbinarize (Multinetwork)" begin + network = WaterModels.parse_file("../test/data/epanet/multinetwork/pump-hw-lps.inp") + network_mn = WaterModels.make_multinetwork(network) + wm = instantiate_model(network_mn, NCWaterModel, build_mn_wf) + unbinarize_mn!(wm, 3) # Transform binary variables to continuous. + @test !JuMP.is_binary(_IM.var(wm, 1, :z_pump, 1)) + @test !JuMP.is_binary(_IM.var(wm, 2, :z_pump, 1)) + @test !JuMP.is_binary(_IM.var(wm, 3, :z_pump, 1)) +end diff --git a/test/wf.jl b/test/wf.jl index f1b2c2c8..74a64c0b 100644 --- a/test/wf.jl +++ b/test/wf.jl @@ -9,25 +9,35 @@ @test result["termination_status"] == LOCALLY_SOLVED @test isapprox(result["solution"]["node"]["1"]["h"], 10.0, rtol=1.0e-3) @test isapprox(result["solution"]["node"]["2"]["p"], 7.89, rtol=1.0e-3) - @test isapprox(result["solution"][name]["1"]["q"], 1.0, rtol=1.0e-3) + @test isapprox(result["solution"]["pipe"]["1"]["q"], 1.0, rtol=1.0e-3) wm = instantiate_model(network, CRDWaterModel, build_wf) result = WaterModels.optimize_model!(wm, optimizer=_make_juniper(wm, ipopt)) @test result["termination_status"] == LOCALLY_SOLVED @test isapprox(result["solution"]["node"]["1"]["h"], 10.0, rtol=1.0e-3) @test result["solution"]["node"]["2"]["p"] <= 7.89 - @test isapprox(result["solution"][name]["1"]["q"], 1.0, rtol=1.0e-3) + @test isapprox(result["solution"]["pipe"]["1"]["q"], 1.0, rtol=1.0e-3) result = run_wf(network, LAWaterModel, cbc, ext=Dict(:pipe_breakpoints=>3)) @test result["termination_status"] == OPTIMAL @test isapprox(result["solution"]["node"]["1"]["h"], 10.0, rtol=1.0e-3) @test isapprox(result["solution"]["node"]["2"]["p"], 7.89, rtol=1.0e-3) - @test isapprox(result["solution"][name]["1"]["q"], 1.0, rtol=1.0e-3) + @test isapprox(result["solution"]["pipe"]["1"]["q"], 1.0, rtol=1.0e-3) result = run_wf(network, LRDWaterModel, cbc, ext=Dict(:pipe_breakpoints=>3)) @test result["termination_status"] == OPTIMAL @test isapprox(result["solution"]["node"]["1"]["h"], 10.0, rtol=1.0e-3) - @test isapprox(result["solution"][name]["1"]["q"], 1.0, rtol=1.0e-3) + @test isapprox(result["solution"]["pipe"]["1"]["q"], 1.0, rtol=1.0e-3) + + result = run_wf(network, QRDWaterModel, cbc, ext=Dict(:pipe_breakpoints=>3)) + @test result["termination_status"] == OPTIMAL + @test isapprox(result["solution"]["node"]["1"]["h"], 10.0, rtol=1.0e-3) + @test isapprox(result["solution"]["pipe"]["1"]["q"], 1.0, rtol=1.0e-3) + + result = run_wf(network, CQRDWaterModel, cbc, ext=Dict(:pipe_breakpoints=>3)) + @test result["termination_status"] == OPTIMAL + @test isapprox(result["solution"]["node"]["1"]["h"], 10.0, rtol=1.0e-3) + @test isapprox(result["solution"]["pipe"]["1"]["q"], 1.0, rtol=1.0e-3) end end @@ -59,9 +69,21 @@ @test result["termination_status"] == OPTIMAL @test isapprox(result["solution"]["node"]["1"]["h"], 10.0, rtol=1.0e-3) @test isapprox(result["solution"]["pipe"]["1"]["q"], 0.9, rtol=1.0e-3) + + wm = instantiate_model(network, QRDWaterModel, build_wf, ext=Dict(:pipe_breakpoints=>3)) + result = WaterModels.optimize_model!(wm, optimizer=_make_juniper(wm, ipopt; register=false)) + @test result["termination_status"] == LOCALLY_SOLVED + @test isapprox(result["solution"]["node"]["1"]["h"], 10.0, rtol=1.0e-3) + @test isapprox(result["solution"]["pipe"]["1"]["q"], 0.9, rtol=1.0e-3) + + wm = instantiate_model(network, CQRDWaterModel, build_wf, ext=Dict(:pipe_breakpoints=>3)) + result = WaterModels.optimize_model!(wm, optimizer=_make_juniper(wm, ipopt; register=false)) + @test result["termination_status"] == LOCALLY_SOLVED + @test isapprox(result["solution"]["node"]["1"]["h"], 10.0, rtol=1.0e-3) + @test isapprox(result["solution"]["pipe"]["1"]["q"], 0.9, rtol=1.0e-3) end - @testset "Head Loss (PRV)" begin + @testset "Head Loss (Regulator)" begin network = WaterModels.parse_file("../test/data/epanet/snapshot/prv-hw-lps.inp") wm = instantiate_model(network, NCWaterModel, build_wf) @@ -95,6 +117,22 @@ @test isapprox(result["solution"]["node"]["3"]["h"], 2.00, rtol=1.0e-3) @test isapprox(result["solution"]["node"]["3"]["p"], 1.00, rtol=1.0e-3) @test isapprox(result["solution"]["pipe"]["1"]["q"], 2.0, rtol=1.0e-3) + + wm = instantiate_model(network, QRDWaterModel, build_wf, ext=Dict(:pipe_breakpoints=>3)) + result = WaterModels.optimize_model!(wm, optimizer=_make_juniper(wm, ipopt; register=false)) + @test result["termination_status"] == LOCALLY_SOLVED + @test isapprox(result["solution"]["node"]["1"]["h"], 10.0, rtol=1.0e-3) + @test isapprox(result["solution"]["node"]["3"]["h"], 2.00, rtol=1.0e-3) + @test isapprox(result["solution"]["node"]["3"]["p"], 1.00, rtol=1.0e-3) + @test isapprox(result["solution"]["pipe"]["1"]["q"], 2.0, rtol=1.0e-3) + + wm = instantiate_model(network, CQRDWaterModel, build_wf, ext=Dict(:pipe_breakpoints=>3)) + result = WaterModels.optimize_model!(wm, optimizer=_make_juniper(wm, ipopt; register=false)) + @test result["termination_status"] == LOCALLY_SOLVED + @test isapprox(result["solution"]["node"]["1"]["h"], 10.0, rtol=1.0e-3) + @test isapprox(result["solution"]["node"]["3"]["h"], 2.00, rtol=1.0e-3) + @test isapprox(result["solution"]["node"]["3"]["p"], 1.00, rtol=1.0e-3) + @test isapprox(result["solution"]["pipe"]["1"]["q"], 2.0, rtol=1.0e-3) end @testset "Head Gain (Pump)" begin @@ -123,6 +161,19 @@ @test result["termination_status"] == OPTIMAL @test isapprox(result["solution"]["node"]["1"]["h"], 10.0, rtol=1.0e-3) @test isapprox(result["solution"]["pump"]["1"]["status"], 1.0, atol=1.0e-3) + + wm = instantiate_model(network, QRDWaterModel, build_wf) + result = WaterModels.optimize_model!(wm, optimizer=_make_juniper(wm, ipopt; register=false)) + @test result["termination_status"] == LOCALLY_SOLVED + @test isapprox(result["solution"]["node"]["1"]["h"], 10.0, rtol=1.0e-3) + @test isapprox(result["solution"]["node"]["2"]["h"], 98.98, rtol=1.0e-3) + @test isapprox(result["solution"]["pump"]["1"]["status"], 1.0, atol=1.0e-3) + + wm = instantiate_model(network, CQRDWaterModel, build_wf) + result = WaterModels.optimize_model!(wm, optimizer=_make_juniper(wm, ipopt; register=false)) + @test result["termination_status"] == LOCALLY_SOLVED + @test isapprox(result["solution"]["node"]["1"]["h"], 10.0, rtol=1.0e-3) + @test isapprox(result["solution"]["pump"]["1"]["status"], 1.0, atol=1.0e-3) end @testset "Hazen-Williams Head Loss (Tank)" begin @@ -152,6 +203,20 @@ @test result["termination_status"] == OPTIMAL @test isapprox(result["solution"]["node"]["1"]["h"], 20.0, rtol=1.0e-3) @test isapprox(result["solution"]["pipe"]["1"]["q"], 1.0, rtol=1.0e-3) + + wm = instantiate_model(network, QRDWaterModel, build_wf, ext=Dict(:pipe_breakpoints=>3)) + result = WaterModels.optimize_model!(wm, optimizer=_make_juniper(wm, ipopt; register=false)) + @test result["termination_status"] == LOCALLY_SOLVED + @test isapprox(result["solution"]["node"]["1"]["h"], 20.0, rtol=1.0e-3) + @test result["solution"]["node"]["2"]["h"] <= 17.89 + @test isapprox(result["solution"]["pipe"]["1"]["q"], 1.0, rtol=1.0e-3) + + wm = instantiate_model(network, CQRDWaterModel, build_wf, ext=Dict(:pipe_breakpoints=>3)) + result = WaterModels.optimize_model!(wm, optimizer=_make_juniper(wm, ipopt; register=false)) + @test result["termination_status"] == LOCALLY_SOLVED + @test isapprox(result["solution"]["node"]["1"]["h"], 20.0, rtol=1.0e-3) + @test result["solution"]["node"]["2"]["h"] <= 17.89 + @test isapprox(result["solution"]["pipe"]["1"]["q"], 1.0, rtol=1.0e-3) end end @@ -169,17 +234,17 @@ end @test isapprox(result["solution"]["nw"]["1"]["node"]["2"]["p"], 7.89, rtol=1.0e-3) @test isapprox(result["solution"]["nw"]["2"]["node"]["2"]["h"], 9.42, rtol=1.0e-3) @test isapprox(result["solution"]["nw"]["3"]["node"]["2"]["h"], 9.84, rtol=1.0e-3) - @test isapprox(result["solution"]["nw"]["1"][name]["1"]["q"], 1.0, rtol=1.0e-3) - @test isapprox(result["solution"]["nw"]["2"][name]["1"]["q"], 0.5, rtol=1.0e-3) - @test isapprox(result["solution"]["nw"]["3"][name]["1"]["q"], 0.25, rtol=1.0e-3) + @test isapprox(result["solution"]["nw"]["1"]["pipe"]["1"]["q"], 1.0, rtol=1.0e-3) + @test isapprox(result["solution"]["nw"]["2"]["pipe"]["1"]["q"], 0.5, rtol=1.0e-3) + @test isapprox(result["solution"]["nw"]["3"]["pipe"]["1"]["q"], 0.25, rtol=1.0e-3) wm = instantiate_model(network, CRDWaterModel, build_mn_wf) result = WaterModels.optimize_model!(wm, optimizer=_make_juniper(wm, ipopt)) @test result["termination_status"] == LOCALLY_SOLVED @test isapprox(result["solution"]["nw"]["2"]["node"]["1"]["h"], 10.0, rtol=1.0e-3) - @test isapprox(result["solution"]["nw"]["1"][name]["1"]["q"], 1.0, rtol=1.0e-3) - @test isapprox(result["solution"]["nw"]["2"][name]["1"]["q"], 0.5, rtol=1.0e-3) - @test isapprox(result["solution"]["nw"]["3"][name]["1"]["q"], 0.25, rtol=1.0e-3) + @test isapprox(result["solution"]["nw"]["1"]["pipe"]["1"]["q"], 1.0, rtol=1.0e-3) + @test isapprox(result["solution"]["nw"]["2"]["pipe"]["1"]["q"], 0.5, rtol=1.0e-3) + @test isapprox(result["solution"]["nw"]["3"]["pipe"]["1"]["q"], 0.25, rtol=1.0e-3) result = run_mn_wf(network, LAWaterModel, cbc, ext=Dict(:pipe_breakpoints=>3)) @test result["termination_status"] == OPTIMAL @@ -187,16 +252,32 @@ end @test isapprox(result["solution"]["nw"]["1"]["node"]["2"]["p"], 7.89, rtol=1.0e-3) @test isapprox(result["solution"]["nw"]["2"]["node"]["2"]["h"], 9.42, rtol=1.0e-3) @test isapprox(result["solution"]["nw"]["3"]["node"]["2"]["h"], 9.84, rtol=1.0e-3) - @test isapprox(result["solution"]["nw"]["1"][name]["1"]["q"], 1.0, rtol=1.0e-3) - @test isapprox(result["solution"]["nw"]["2"][name]["1"]["q"], 0.5, rtol=1.0e-3) - @test isapprox(result["solution"]["nw"]["3"][name]["1"]["q"], 0.25, rtol=1.0e-3) + @test isapprox(result["solution"]["nw"]["1"]["pipe"]["1"]["q"], 1.0, rtol=1.0e-3) + @test isapprox(result["solution"]["nw"]["2"]["pipe"]["1"]["q"], 0.5, rtol=1.0e-3) + @test isapprox(result["solution"]["nw"]["3"]["pipe"]["1"]["q"], 0.25, rtol=1.0e-3) result = run_mn_wf(network, LRDWaterModel, cbc, ext=Dict(:pipe_breakpoints=>3)) @test result["termination_status"] == OPTIMAL @test isapprox(result["solution"]["nw"]["2"]["node"]["1"]["h"], 10.0, rtol=1.0e-3) - @test isapprox(result["solution"]["nw"]["1"][name]["1"]["q"], 1.0, rtol=1.0e-3) - @test isapprox(result["solution"]["nw"]["2"][name]["1"]["q"], 0.5, rtol=1.0e-3) - @test isapprox(result["solution"]["nw"]["3"][name]["1"]["q"], 0.25, rtol=1.0e-3) + @test isapprox(result["solution"]["nw"]["1"]["pipe"]["1"]["q"], 1.0, rtol=1.0e-3) + @test isapprox(result["solution"]["nw"]["2"]["pipe"]["1"]["q"], 0.5, rtol=1.0e-3) + @test isapprox(result["solution"]["nw"]["3"]["pipe"]["1"]["q"], 0.25, rtol=1.0e-3) + + wm = instantiate_model(network, QRDWaterModel, build_mn_wf, ext=Dict(:pipe_breakpoints=>3)) + result = WaterModels.optimize_model!(wm, optimizer=_make_juniper(wm, ipopt; register=false)) + @test result["termination_status"] == LOCALLY_SOLVED + @test isapprox(result["solution"]["nw"]["2"]["node"]["1"]["h"], 10.0, rtol=1.0e-3) + @test isapprox(result["solution"]["nw"]["1"]["pipe"]["1"]["q"], 1.0, rtol=1.0e-3) + @test isapprox(result["solution"]["nw"]["2"]["pipe"]["1"]["q"], 0.5, rtol=1.0e-3) + @test isapprox(result["solution"]["nw"]["3"]["pipe"]["1"]["q"], 0.25, rtol=1.0e-3) + + wm = instantiate_model(network, CQRDWaterModel, build_mn_wf, ext=Dict(:pipe_breakpoints=>3)) + result = WaterModels.optimize_model!(wm, optimizer=_make_juniper(wm, ipopt; register=false)) + @test result["termination_status"] == LOCALLY_SOLVED + @test isapprox(result["solution"]["nw"]["2"]["node"]["1"]["h"], 10.0, rtol=1.0e-3) + @test isapprox(result["solution"]["nw"]["1"]["pipe"]["1"]["q"], 1.0, rtol=1.0e-3) + @test isapprox(result["solution"]["nw"]["2"]["pipe"]["1"]["q"], 0.5, rtol=1.0e-3) + @test isapprox(result["solution"]["nw"]["3"]["pipe"]["1"]["q"], 0.25, rtol=1.0e-3) end end @@ -229,9 +310,21 @@ end @test result["termination_status"] == OPTIMAL @test isapprox(result["solution"]["nw"]["1"]["node"]["1"]["h"], 10.0, rtol=1.0e-3) @test isapprox(result["solution"]["nw"]["2"]["pipe"]["1"]["q"], 0.9, rtol=1.0e-3) + + wm = instantiate_model(network, QRDWaterModel, build_mn_wf, ext=Dict(:pipe_breakpoints=>3)) + result = WaterModels.optimize_model!(wm, optimizer=_make_juniper(wm, ipopt; register=false)) + @test result["termination_status"] == LOCALLY_SOLVED + @test isapprox(result["solution"]["nw"]["1"]["node"]["1"]["h"], 10.0, rtol=1.0e-3) + @test isapprox(result["solution"]["nw"]["3"]["pipe"]["1"]["q"], 0.9, rtol=1.0e-3) + + wm = instantiate_model(network, CQRDWaterModel, build_mn_wf, ext=Dict(:pipe_breakpoints=>3)) + result = WaterModels.optimize_model!(wm, optimizer=_make_juniper(wm, ipopt; register=false)) + @test result["termination_status"] == LOCALLY_SOLVED + @test isapprox(result["solution"]["nw"]["1"]["node"]["1"]["h"], 10.0, rtol=1.0e-3) + @test isapprox(result["solution"]["nw"]["3"]["pipe"]["1"]["q"], 0.9, rtol=1.0e-3) end - @testset "Head Loss (PRV)" begin + @testset "Head Loss (Regulator)" begin network = WaterModels.parse_file("../test/data/epanet/multinetwork/prv-hw-lps.inp") network = WaterModels.make_multinetwork(network) @@ -245,8 +338,8 @@ end @test isapprox(result["solution"]["nw"]["2"]["node"]["1"]["h"], 10.0, rtol=1.0e-3) @test isapprox(result["solution"]["nw"]["1"]["pipe"]["1"]["q"], 2.00, rtol=1.0e-3) @test isapprox(result["solution"]["nw"]["3"]["pipe"]["1"]["q"], 0.50, rtol=1.0e-3) - @test isapprox(result["solution"]["nw"]["1"]["pressure_reducing_valve"]["2"]["q"], 1.00, rtol=1.0e-3) - @test isapprox(result["solution"]["nw"]["3"]["pressure_reducing_valve"]["2"]["q"], 0.25, rtol=1.0e-3) + @test isapprox(result["solution"]["nw"]["1"]["regulator"]["2"]["q"], 1.00, rtol=1.0e-3) + @test isapprox(result["solution"]["nw"]["3"]["regulator"]["2"]["q"], 0.25, rtol=1.0e-3) wm = instantiate_model(network, CRDWaterModel, build_mn_wf) result = WaterModels.optimize_model!(wm, optimizer=_make_juniper(wm, ipopt)) @@ -258,8 +351,8 @@ end @test isapprox(result["solution"]["nw"]["2"]["node"]["1"]["h"], 10.0, rtol=1.0e-3) @test isapprox(result["solution"]["nw"]["1"]["pipe"]["1"]["q"], 2.00, rtol=1.0e-3) @test isapprox(result["solution"]["nw"]["3"]["pipe"]["1"]["q"], 0.50, rtol=1.0e-3) - @test isapprox(result["solution"]["nw"]["1"]["pressure_reducing_valve"]["2"]["q"], 1.00, rtol=1.0e-3) - @test isapprox(result["solution"]["nw"]["3"]["pressure_reducing_valve"]["2"]["q"], 0.25, rtol=1.0e-3) + @test isapprox(result["solution"]["nw"]["1"]["regulator"]["2"]["q"], 1.00, rtol=1.0e-3) + @test isapprox(result["solution"]["nw"]["3"]["regulator"]["2"]["q"], 0.25, rtol=1.0e-3) result = run_mn_wf(network, LAWaterModel, cbc, ext=Dict(:pipe_breakpoints=>3)) @test result["termination_status"] == OPTIMAL @@ -270,20 +363,35 @@ end @test isapprox(result["solution"]["nw"]["2"]["node"]["1"]["h"], 10.0, rtol=1.0e-3) @test isapprox(result["solution"]["nw"]["1"]["pipe"]["1"]["q"], 2.00, rtol=1.0e-3) @test isapprox(result["solution"]["nw"]["3"]["pipe"]["1"]["q"], 0.50, rtol=1.0e-3) - @test isapprox(result["solution"]["nw"]["1"]["pressure_reducing_valve"]["2"]["q"], 1.00, rtol=1.0e-3) - @test isapprox(result["solution"]["nw"]["3"]["pressure_reducing_valve"]["2"]["q"], 0.25, rtol=1.0e-3) + @test isapprox(result["solution"]["nw"]["1"]["regulator"]["2"]["q"], 1.00, rtol=1.0e-3) + @test isapprox(result["solution"]["nw"]["3"]["regulator"]["2"]["q"], 0.25, rtol=1.0e-3) result = run_mn_wf(network, LRDWaterModel, cbc, ext=Dict(:pipe_breakpoints=>3)) @test result["termination_status"] == OPTIMAL - @test result["solution"]["nw"]["1"]["node"]["2"]["h"] <= 2.39 - @test result["solution"]["nw"]["2"]["node"]["2"]["h"] <= 7.90 - @test result["solution"]["nw"]["3"]["node"]["2"]["h"] <= 9.42 @test isapprox(result["solution"]["nw"]["2"]["node"]["3"]["h"], 2.00, rtol=1.0e-3) @test isapprox(result["solution"]["nw"]["2"]["node"]["1"]["h"], 10.0, rtol=1.0e-3) @test isapprox(result["solution"]["nw"]["1"]["pipe"]["1"]["q"], 2.00, rtol=1.0e-3) @test isapprox(result["solution"]["nw"]["3"]["pipe"]["1"]["q"], 0.50, rtol=1.0e-3) - @test isapprox(result["solution"]["nw"]["1"]["pressure_reducing_valve"]["2"]["q"], 1.00, rtol=1.0e-3) - @test isapprox(result["solution"]["nw"]["3"]["pressure_reducing_valve"]["2"]["q"], 0.25, rtol=1.0e-3) + @test isapprox(result["solution"]["nw"]["1"]["regulator"]["2"]["q"], 1.00, rtol=1.0e-3) + @test isapprox(result["solution"]["nw"]["3"]["regulator"]["2"]["q"], 0.25, rtol=1.0e-3) + + result = run_mn_wf(network, QRDWaterModel, cbc, ext=Dict(:pipe_breakpoints=>3)) + @test result["termination_status"] == OPTIMAL + @test isapprox(result["solution"]["nw"]["2"]["node"]["3"]["h"], 2.00, rtol=1.0e-3) + @test isapprox(result["solution"]["nw"]["2"]["node"]["1"]["h"], 10.0, rtol=1.0e-3) + @test isapprox(result["solution"]["nw"]["1"]["pipe"]["1"]["q"], 2.00, rtol=1.0e-3) + @test isapprox(result["solution"]["nw"]["3"]["pipe"]["1"]["q"], 0.50, rtol=1.0e-3) + @test isapprox(result["solution"]["nw"]["1"]["regulator"]["2"]["q"], 1.00, rtol=1.0e-3) + @test isapprox(result["solution"]["nw"]["3"]["regulator"]["2"]["q"], 0.25, rtol=1.0e-3) + + result = run_mn_wf(network, CQRDWaterModel, cbc, ext=Dict(:pipe_breakpoints=>3)) + @test result["termination_status"] == OPTIMAL + @test isapprox(result["solution"]["nw"]["2"]["node"]["3"]["h"], 2.00, rtol=1.0e-3) + @test isapprox(result["solution"]["nw"]["2"]["node"]["1"]["h"], 10.0, rtol=1.0e-3) + @test isapprox(result["solution"]["nw"]["1"]["pipe"]["1"]["q"], 2.00, rtol=1.0e-3) + @test isapprox(result["solution"]["nw"]["3"]["pipe"]["1"]["q"], 0.50, rtol=1.0e-3) + @test isapprox(result["solution"]["nw"]["1"]["regulator"]["2"]["q"], 1.00, rtol=1.0e-3) + @test isapprox(result["solution"]["nw"]["3"]["regulator"]["2"]["q"], 0.25, rtol=1.0e-3) end @testset "Head Gain (Pump)" begin @@ -320,6 +428,23 @@ end @test isapprox(result["solution"]["nw"]["1"]["pump"]["1"]["q"], 0.125, rtol=1.0e-3) @test isapprox(result["solution"]["nw"]["3"]["pump"]["1"]["q"], 0.03125, rtol=1.0e-3) @test isapprox(result["solution"]["nw"]["1"]["pump"]["1"]["status"], 1.0, atol=1.0e-3) + + wm = instantiate_model(network, QRDWaterModel, build_mn_wf) + result = WaterModels.optimize_model!(wm, optimizer=_make_juniper(wm, ipopt; register=false)) + @test result["termination_status"] == LOCALLY_SOLVED + @test isapprox(result["solution"]["nw"]["1"]["pump"]["1"]["q"], 0.125, rtol=1.0e-3) + @test isapprox(result["solution"]["nw"]["3"]["pump"]["1"]["q"], 0.03125, rtol=1.0e-3) + @test isapprox(result["solution"]["nw"]["1"]["pump"]["1"]["status"], 1.0, atol=1.0e-3) + @test isapprox(result["solution"]["nw"]["1"]["pump"]["1"]["g"], 88.98, rtol=1.0e-3) + @test isapprox(result["solution"]["nw"]["3"]["pump"]["1"]["g"], 99.59, rtol=1.0e-2) + + wm = instantiate_model(network, CQRDWaterModel, build_mn_wf) + result = WaterModels.optimize_model!(wm, optimizer=_make_juniper(wm, ipopt; register=false)) + @test result["termination_status"] == LOCALLY_SOLVED + @test isapprox(result["solution"]["nw"]["1"]["pump"]["1"]["q"], 0.125, rtol=1.0e-3) + @test isapprox(result["solution"]["nw"]["3"]["pump"]["1"]["q"], 0.03125, rtol=1.0e-3) + @test isapprox(result["solution"]["nw"]["1"]["pump"]["1"]["status"], 1.0, atol=1.0e-3) + @test result["solution"]["nw"]["1"]["pump"]["1"]["g"] <= 88.99 end @testset "Hazen-Williams Head Loss (Tank)" begin @@ -361,5 +486,25 @@ end @test isapprox(result["solution"]["nw"]["3"]["node"]["1"]["h"], 25.62, rtol=1.0e-3) @test isapprox(result["solution"]["nw"]["1"]["pipe"]["1"]["q"], 0.500, rtol=1.0e-3) @test isapprox(result["solution"]["nw"]["3"]["pipe"]["1"]["q"], 0.125, rtol=1.0e-3) + + wm = instantiate_model(network, QRDWaterModel, build_mn_wf, ext=Dict(:pipe_breakpoints=>4)) + result = WaterModels.optimize_model!(wm, optimizer=_make_juniper(wm, ipopt; register=false)) + @test result["termination_status"] == LOCALLY_SOLVED + @test isapprox(result["solution"]["nw"]["1"]["node"]["1"]["h"], 60.00, rtol=1.0e-3) + @test result["solution"]["nw"]["1"]["node"]["2"]["h"] <= 59.42 + @test isapprox(result["solution"]["nw"]["3"]["node"]["1"]["h"], 25.62, rtol=1.0e-3) + @test result["solution"]["nw"]["3"]["node"]["2"]["h"] <= 25.58 + @test isapprox(result["solution"]["nw"]["1"]["pipe"]["1"]["q"], 0.500, rtol=1.0e-3) + @test isapprox(result["solution"]["nw"]["3"]["pipe"]["1"]["q"], 0.125, rtol=1.0e-3) + + wm = instantiate_model(network, CQRDWaterModel, build_mn_wf, ext=Dict(:pipe_breakpoints=>4)) + result = WaterModels.optimize_model!(wm, optimizer=_make_juniper(wm, ipopt; register=false)) + @test result["termination_status"] == LOCALLY_SOLVED + @test isapprox(result["solution"]["nw"]["1"]["node"]["1"]["h"], 60.00, rtol=1.0e-3) + @test result["solution"]["nw"]["1"]["node"]["2"]["h"] <= 59.42 + @test isapprox(result["solution"]["nw"]["3"]["node"]["1"]["h"], 25.62, rtol=1.0e-3) + @test result["solution"]["nw"]["3"]["node"]["2"]["h"] <= 25.58 + @test isapprox(result["solution"]["nw"]["1"]["pipe"]["1"]["q"], 0.500, rtol=1.0e-3) + @test isapprox(result["solution"]["nw"]["3"]["pipe"]["1"]["q"], 0.125, rtol=1.0e-3) end end