Merge pull request #103 from lanl-ansi/obbt

ADD: Optimization-based bound tightening utility

.travis.yml

@@ -4,15 +4,15 @@ os:
   - osx
 julia:
   - 1.0
-  - 1.5
+  - 1
   - nightly
 codecov: true
 jobs:
   allow_failures:
     - julia: nightly
   include:
     - stage: "Documentation"
-      julia: 1.5
+      julia: 1
       os: linux
       script:
         - julia --project=docs/ -e 'using Pkg; Pkg.instantiate(); Pkg.develop(PackageSpec(path=pwd()))'

CHANGELOG.md

@@ -9,6 +9,7 @@ WaterModels.jl Change Log
 - Fixed various bugs related to component indices.
 - Added new constraints to tighten `owf` formulations.
 - Added utility function for "unbinarizing" indicator variables.
+- Added utility function for optimization-based bound tightening.
 - Simplified integration tests and removed previous tests.
 
 ### v0.3.0

src/WaterModels.jl

@@ -67,6 +67,7 @@ include("prob/owf.jl")
 include("prob/des.jl")
 
 include("util/unbinarize.jl")
+include("util/obbt.jl")
 
 include("core/export.jl")
 end

src/core/base.jl

@@ -154,6 +154,15 @@ 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_reservoir] = filter(x -> x.second["dispatchable"], ref[:reservoir])
+        ref[:nondispatchable_reservoir] = filter(x -> !x.second["dispatchable"], ref[:reservoir])
+        ref[:dispatchable_tank] = filter(x -> x.second["dispatchable"], ref[:tank])
+        ref[:nondispatchable_tank] = filter(x -> !x.second["dispatchable"], ref[:tank])
+
+        # 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])

src/core/constraint.jl

@@ -25,7 +25,8 @@ function constraint_flow_conservation(
     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}, demand::Float64)
+    reservoirs::Array{Int64,1}, tanks::Array{Int64,1},
+    dispatchable_junctions::Array{Int64,1}, fixed_demand::Float64)
     # Collect flow variable references per component.
     q_check_valve = var(wm, n, :q_check_valve)
     q_pipe, q_pump = var(wm, n, :q_pipe), var(wm, n, :q_pump)
@@ -43,7 +44,8 @@ function constraint_flow_conservation(
          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) + demand)
+         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
 end

src/core/constraint_template.jl

@@ -45,8 +45,12 @@ function constraint_flow_conservation(wm::AbstractWaterModel, i::Int; nw::Int=wm
     junctions = ref(wm, nw, :node_junction, i) # Junctions attached to node `i`.
 
     # Sum the constant demands required at node `i`.
-    demands = [ref(wm, nw, :junction, k)["demand"] for k in junctions]
-    net_demand = length(demands) > 0 ? sum(demands) : 0.0
+    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]
+    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)
 
     # Initialize the flow conservation constraint dictionary entry.
     _initialize_con_dict(wm, :flow_conservation, nw=nw)
@@ -55,7 +59,48 @@ function constraint_flow_conservation(wm::AbstractWaterModel, i::Int; nw::Int=wm
     constraint_flow_conservation(
         wm, nw, i, check_valve_fr, check_valve_to, pipe_fr, pipe_to, pump_fr, pump_to,
         pressure_reducing_valve_fr, pressure_reducing_valve_to, shutoff_valve_fr,
-        shutoff_valve_to, reservoirs, tanks, net_demand)
+        shutoff_valve_to, reservoirs, tanks, dispatchable_junctions, 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)
+    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)
+
+    # 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)
+
+    # Get the in degree of node `i`.
+    in_length = length(check_valve_to) + length(pipe_to) + length(pump_to) +
+                length(pressure_reducing_valve_to) + length(shutoff_valve_to)
+
+    # Get the out degree of node `i`.
+    out_length = length(check_valve_fr) + length(pipe_fr) + length(pump_fr) +
+                 length(pressure_reducing_valve_fr) + length(shutoff_valve_fr)
+
+    # 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, pump_fr, pump_to,
+            pressure_reducing_valve_fr, pressure_reducing_valve_to, shutoff_valve_fr,
+            shutoff_valve_to)
+    end
 end
 
 
@@ -86,10 +131,13 @@ end
 
 ### Reservoir Constraints ###
 function constraint_reservoir_head(wm::AbstractWaterModel, i::Int; nw::Int=wm.cnw)
-    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)
+    # 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)
+    end
 end
 
 
@@ -107,7 +155,7 @@ function constraint_source_directionality(wm::AbstractWaterModel, i::Int; nw::In
     shutoff_valve_to = _collect_comps_to(wm, i, :shutoff_valve; nw=nw)
 
     # Initialize the source flow constraint dictionary entry.
-    _initialize_con_dict(wm, :source_flow, nw=nw)
+    _initialize_con_dict(wm, :source_directionality, nw=nw)
 
     # Add the source flow directionality constraint.
     constraint_source_directionality(
@@ -123,13 +171,15 @@ function constraint_tank_state(wm::AbstractWaterModel, i::Int; nw::Int=wm.cnw)
         Memento.error(_LOGGER, "Tank states cannot be controlled without a time step.")
     end
 
-    tank = ref(wm, nw, :tank, i)
-    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)
+    # Only set the tank state if the tank is nondispatchable.
+    if !ref(wm, nw, :tank, i)["dispatchable"]
+        tank = ref(wm, nw, :tank, i)
+        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)
+    end
 end
 
 
@@ -138,10 +188,13 @@ function constraint_tank_state(wm::AbstractWaterModel, i::Int, nw_1::Int, nw_2::
         Memento.error(_LOGGER, "Tank states cannot be controlled without a time step.")
     end
 
-    # 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))
+    # Only set the tank state if the tank is nondispatchable.
+    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))
+    end
 end
 
 

src/core/data.jl

@@ -115,6 +115,18 @@ function calc_resistance_costs_hw(pipes::Dict{Int, <:Any})
 end
 
 
+function make_tank_start_dispatchable!(data::Dict{String,<:Any})
+    if _IM.ismultinetwork(data)
+        nw_ids = sort(collect(keys(data["nw"])))
+        start_nw = string(sort([parse(Int, i) for i in nw_ids])[1])
+
+        for (i, tank) in data["nw"][start_nw]["tank"]
+            tank["dispatchable"] = true
+        end
+    end
+end
+
+
 function calc_resistance_costs_dw(pipes::Dict{Int, <:Any}, viscosity::Float64)
     # Create placeholder costs dictionary.
     costs = Dict([(a, Array{Float64, 1}()) for a in keys(pipes)])

src/core/objective.jl

@@ -12,15 +12,6 @@ function objective_wf(wm::AbstractWaterModel)
     JuMP.set_objective_sense(wm.model, _MOI.FEASIBILITY_SENSE)
 end
 
-"""
-    objective_owf(wm::AbstractWaterModel)
-
-Sets the objective function for optimal water flow (owf) problem
-specifications. By default, only feasibility must be satisfied.
-"""
-function objective_owf(wm::AbstractWaterModel) 
-    JuMP.set_objective_sense(wm.model, _MOI.FEASIBILITY_SENSE)
-end
 
 """
     objective_des(wm::AbstractWaterModel)

src/core/ref.jl

@@ -31,6 +31,20 @@ function _get_function_from_head_curve(head_curve::Array{Tuple{Float64,Float64}}
 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))
+
+    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"])
+    end
+
+    return demand_min, demand_max
+end
+
+
 function calc_head_bounds(wm::AbstractWaterModel, n::Int=wm.cnw)
     # Compute the maximum elevation of all nodes in the network.
     max_head = maximum(node["elevation"] for (i, node) in ref(wm, n, :node))
@@ -52,16 +66,6 @@ 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)
-        # Override the minimum head value at node i, if additional data is present.
-        if haskey(node, "minimumHead")
-            head_min[i] = max(head_min[i], node["minimumHead"])
-        end
-
-        # Override the maximum head value at node i, if additional data is present.
-        if haskey(node, "maximumHead")
-            head_max[i] = min(head_max[i], node["maximumHead"])
-        end
-
         num_junctions = length(ref(wm, n, :node_junction, i))
         num_reservoirs = length(ref(wm, n, :node_reservoir, i))
         num_tanks = length(ref(wm, n, :node_tank, i))
@@ -75,8 +79,14 @@ function calc_head_bounds(wm::AbstractWaterModel, n::Int=wm.cnw)
     for (i, reservoir) in ref(wm, n, :reservoir)
         # Fix head values at nodes with reservoirs to predefined values.
         node_id = reservoir["node"]
-        head_min[node_id] = ref(wm, n, :node, node_id)["head"]
-        head_max[node_id] = ref(wm, n, :node, node_id)["head"]
+
+        if !reservoir["dispatchable"]
+            head_min[node_id] = ref(wm, n, :node, node_id)["head"]
+            head_max[node_id] = ref(wm, n, :node, node_id)["head"]
+        else
+            head_min[node_id] = ref(wm, n, :node, node_id)["h_min"]
+            head_max[node_id] = ref(wm, n, :node, node_id)["h_max"]
+        end
     end
 
     for (i, tank) in ref(wm, n, :tank)
@@ -93,6 +103,11 @@ function calc_head_bounds(wm::AbstractWaterModel, n::Int=wm.cnw)
         head_min[node_to] = min(head_min[node_to], h_setting)
     end
 
+    for (i, node) in ref(wm, n, :node)
+        haskey(node, "h_min") && (head_min[i] = max(head_min[i], node["h_min"]))
+        haskey(node, "h_max") && (head_max[i] = min(head_max[i], node["h_max"]))
+    end
+
     # Return the dictionaries of lower and upper bounds.
     return head_min, head_max
 end
@@ -102,8 +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)
 
-    junctions = values(ref(wm, n, :junction))
-    sum_demand = sum(abs(junction["demand"]) for junction in junctions)
+    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
 
     if :time_step in keys(ref(wm, n))
         time_step = ref(wm, n, :time_step)
@@ -144,15 +161,20 @@ function calc_flow_bounds(wm::AbstractWaterModel, n::Int=wm.cnw)
                     ub[name][a][r_id] = min(ub[name][a][r_id], rate_bound)
                 end
             end
+
+            haskey(comp, "q_min") && (lb[name][a][1] = max(lb[name][a][1], comp["q_min"]))
+            haskey(comp, "q_max") && (ub[name][a][1] = min(ub[name][a][1], comp["q_max"]))
         end
     end
 
     lb["pressure_reducing_valve"] = Dict{Int,Array{Float64}}()
     ub["pressure_reducing_valve"] = Dict{Int,Array{Float64}}()
 
     for (a, prv) in ref(wm, n, :pressure_reducing_valve)
-        lb["pressure_reducing_valve"][a] = [0.0]
-        ub["pressure_reducing_valve"][a] = [sum_demand]
+        name = "pressure_reducing_valve"
+        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"]))
     end
 
     lb["pump"], ub["pump"] = Dict{Int,Array{Float64}}(), Dict{Int,Array{Float64}}()
@@ -162,7 +184,9 @@ function calc_flow_bounds(wm::AbstractWaterModel, n::Int=wm.cnw)
         c = _get_function_from_head_curve(head_curve)
         q_max = (-c[2] + sqrt(c[2]^2 - 4.0*c[1]*c[3])) * inv(2.0*c[1])
         q_max = max(q_max, (-c[2] - sqrt(c[2]^2 - 4.0*c[1]*c[3])) * inv(2.0*c[1]))
-        lb["pump"][a], ub["pump"][a] = [[0.0], [q_max]]
+        lb["pump"][a], ub["pump"][a] = [[0.0], [min(sum_demand, q_max)]]
+        haskey(pump, "q_min") && (lb["pump"][a][1] = max(lb["pump"][a][1], pump["q_min"]))
+        haskey(pump, "q_max") && (ub["pump"][a][1] = min(ub["pump"][a][1], pump["q_max"]))
     end
 
     return lb, ub

src/core/variable.jl

@@ -68,6 +68,25 @@ function variable_reservoir(wm::AbstractWaterModel; nw::Int=wm.cnw, report::Bool
     report && sol_component_value(wm, nw, :reservoir, :qr, ids(wm, nw, :reservoir), qr)
 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"))
+
+    if bounded
+        demand_lb, 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])
+        end
+    end
+
+    report && sol_component_value(wm, nw, :junction, :demand, ids(wm, nw, :dispatchable_junction), 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)
@@ -77,6 +96,7 @@ function variable_tank(wm::AbstractWaterModel; nw::Int=wm.cnw, report::Bool=true
 
     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"]) -