diff --git a/.travis.yml b/.travis.yml
index ed6b4084..8d242f25 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -4,7 +4,7 @@ os:
   - osx
 julia:
   - 1.0
-  - 1.5
+  - 1
   - nightly
 codecov: true
 jobs:
@@ -12,7 +12,7 @@ jobs:
     - 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()))'
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 02cb8e58..56ebdc6d 100644
--- a/CHANGELOG.md
+++ b/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
diff --git a/src/WaterModels.jl b/src/WaterModels.jl
index 53d41f15..d912ed83 100644
--- a/src/WaterModels.jl
+++ b/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
diff --git a/src/core/base.jl b/src/core/base.jl
index 8eaa0a43..350a1972 100644
--- a/src/core/base.jl
+++ b/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])
diff --git a/src/core/constraint.jl b/src/core/constraint.jl
index e79e9d1f..27cef9b4 100644
--- a/src/core/constraint.jl
+++ b/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
diff --git a/src/core/constraint_template.jl b/src/core/constraint_template.jl
index 5b0c7e95..3e44f060 100644
--- a/src/core/constraint_template.jl
+++ b/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
 
 
diff --git a/src/core/data.jl b/src/core/data.jl
index 888afc40..a0b977fc 100644
--- a/src/core/data.jl
+++ b/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)])
diff --git a/src/core/objective.jl b/src/core/objective.jl
index ae612a98..d857c55f 100644
--- a/src/core/objective.jl
+++ b/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)
diff --git a/src/core/ref.jl b/src/core/ref.jl
index 70cb026c..ef53ae8c 100644
--- a/src/core/ref.jl
+++ b/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,6 +161,9 @@ 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
 
@@ -151,8 +171,10 @@ function calc_flow_bounds(wm::AbstractWaterModel, n::Int=wm.cnw)
     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
diff --git a/src/core/variable.jl b/src/core/variable.jl
index 50a16e3f..24223da7 100644
--- a/src/core/variable.jl
+++ b/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"]) -
diff --git a/src/form/directed_flow.jl b/src/form/directed_flow.jl
index 4ceceb2d..ae0da6bf 100644
--- a/src/form/directed_flow.jl
+++ b/src/form/directed_flow.jl
@@ -190,12 +190,8 @@ function constraint_check_valve_common(wm::AbstractDirectedModel, n::Int, a::Int
     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)
 
-    # Ensure that negative flow is fixed to zero.
-    JuMP.fix(qn, 0.0; force=true)
-
-    # If the check valve is open, flow must be appreciably nonnegative.
+    # If the check valve is closed, flow must be zero.
     c_1 = JuMP.@constraint(wm.model, qp <= JuMP.upper_bound(qp) * z)
-    c_2 = JuMP.@constraint(wm.model, qp >= _q_eps * z)
 
     # Get common head variables and associated data.
     h_i, h_j = var(wm, n, :h, node_fr), var(wm, n, :h, node_to)
@@ -203,23 +199,23 @@ function constraint_check_valve_common(wm::AbstractDirectedModel, n::Int, a::Int
     dhp_ub, dhn_ub = JuMP.upper_bound(dhp), JuMP.upper_bound(dhn)
 
     # When the check valve is closed, positive head loss is not possible.
-    c_3 = JuMP.@constraint(wm.model, dhp <= dhp_ub * z)
+    c_2 = JuMP.@constraint(wm.model, dhp <= dhp_ub * z)
 
     # When the check valve is open, negative head loss is not possible.
-    c_4 = JuMP.@constraint(wm.model, dhn <= dhn_ub * (1.0 - z))
+    c_3 = JuMP.@constraint(wm.model, dhn <= dhn_ub * (1.0 - z))
 
     # Constrain head differences based on direction.
-    c_5 = JuMP.@constraint(wm.model, dhp <= dhp_ub * y)
-    c_6 = JuMP.@constraint(wm.model, dhn <= dhn_ub * (1.0 - y))
+    c_4 = JuMP.@constraint(wm.model, dhp <= dhp_ub * y)
+    c_5 = JuMP.@constraint(wm.model, dhn <= dhn_ub * (1.0 - y))
 
     # Constrain directed flows based on direction.
-    c_7 = JuMP.@constraint(wm.model, qp <= JuMP.upper_bound(qp) * y)
+    c_6 = JuMP.@constraint(wm.model, qp <= JuMP.upper_bound(qp) * y)
 
     # Relate head differences with head variables
-    c_8 = JuMP.@constraint(wm.model, dhp - dhn == h_i - h_j)
+    c_7 = JuMP.@constraint(wm.model, dhp - dhn == h_i - h_j)
 
     # 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, c_8])
+    append!(con(wm, n, :check_valve, a), [c_1, c_2, c_3, c_4, c_5, c_6, c_7])
 end
 
 
@@ -286,9 +282,6 @@ function constraint_prv_common(wm::AbstractDirectedModel, n::Int, a::Int, node_f
     y = var(wm, n, :y_pressure_reducing_valve, a)
     z = var(wm, n, :z_pressure_reducing_valve, a)
 
-    # Ensure that negative flow is fixed to zero.
-    JuMP.fix(qn, 0.0; force=true)
-
     # 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)
@@ -326,7 +319,7 @@ function constraint_pump_head_gain_lb(wm::AbstractDirectedModel, n::Int, a::Int,
     lambda = var(wm, n, :lambda_pump)
 
     # Add a constraint that lower-bounds the head gain variable.
-    breakpoints = range(_q_eps, stop=JuMP.upper_bound(qp), length=pump_breakpoints)
+    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)
@@ -340,9 +333,6 @@ function constraint_pump_common(wm::AbstractDirectedModel, n::Int, a::Int, node_
     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)
 
-    # Ensure that negative flow is fixed to zero.
-    JuMP.fix(var(wm, n, :qn_pump, a), 0.0; force=true)
-
     # 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)
@@ -456,6 +446,48 @@ function constraint_resistance_selection_des(wm::AbstractDirectedModel, n::Int,
 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},
+    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})
+    # Collect direction variable references per component.
+    y_check_valve = var(wm, n, :y_check_valve)
+    y_pipe, y_pump = var(wm, n, :y_pipe), 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)
+
+    y_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_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_out_length = length(check_valve_fr) + length(pipe_fr) + length(pump_fr) +
+                   length(pressure_reducing_valve_fr) + length(shutoff_valve_fr)
+
+    y_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_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))
+
+    y_in_length = length(check_valve_to) + length(pipe_to) + length(pump_to) +
+                  length(pressure_reducing_valve_to) + length(shutoff_valve_to)
+
+    # Add the node directionality constraint.
+    if y_out_length == 1 && y_in_length == 1
+        c = JuMP.@constraint(wm.model, y_out - y_in == 0.0)
+        con(wm, n, :node_directionality)[i] = c
+    elseif y_in_length + y_out_length == 2 && y_in_length*y_out_length == 0
+        c = JuMP.@constraint(wm.model, y_out + y_in == 1.0)
+        con(wm, n, :node_directionality)[i] = c
+    end
+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},
@@ -486,7 +518,7 @@ function constraint_source_directionality(
 
     # Add the source flow direction constraint.
     c = JuMP.@constraint(wm.model, y_out - y_in >= 1.0 - y_in_length)
-    con(wm, n, :source_flow)[i] = c
+    con(wm, n, :source_directionality)[i] = c
 end
 
 
diff --git a/src/form/micp.jl b/src/form/micp.jl
index 5b1b560a..c8545681 100644
--- a/src/form/micp.jl
+++ b/src/form/micp.jl
@@ -32,7 +32,7 @@ function constraint_pump_head_gain(wm::AbstractMICPModel, n::Int, a::Int, node_f
 
     # 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_expr >= g) # Concavified.
+    c = JuMP.@constraint(wm.model, g <= g_expr) # Concavified.
 
     # Append the constraint array.
     append!(con(wm, n, :head_gain, a), [c])
@@ -52,7 +52,7 @@ function constraint_pump_head_gain(wm::AbstractMICPModel, n::Int, a::Int, node_f
 
     # Add a constraint for the flow piecewise approximation.
     qp_ub = JuMP.upper_bound(qp)
-    breakpoints = range(_q_eps, stop=qp_ub, length=pump_breakpoints)
+    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)
 
@@ -80,10 +80,9 @@ function constraint_check_valve_head_loss(wm::AbstractMICPModel, n::Int, a::Int,
     # 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))
-    c_n = JuMP.@NLconstraint(wm.model, dhn <= inv(L) * dhn_ub * (1.0 - z))
 
     # Append the constraint array.
-    append!(con(wm, n, :head_loss)[a], [c_p, c_n])
+    append!(con(wm, n, :head_loss)[a], [c_p])
 end
 
 function constraint_shutoff_valve_head_loss(wm::AbstractMICPModel, n::Int, a::Int, node_fr::Int, node_to::Int, L::Float64, r::Float64)
@@ -154,7 +153,7 @@ function objective_owf(wm::AbstractMICPModel)
                 qp_ub = JuMP.upper_bound(qp)
 
                 # Generate a set of uniform flow and cubic function breakpoints.
-                breakpoints = range(_q_eps, stop=qp_ub, length=pump_breakpoints)
+                breakpoints = range(0.0, stop=qp_ub, length=pump_breakpoints)
                 f = _calc_cubic_flow_values(collect(breakpoints), curve_fun)
 
                 # Get pump efficiency data.
diff --git a/src/form/milp.jl b/src/form/milp.jl
index e37bf874..69710642 100644
--- a/src/form/milp.jl
+++ b/src/form/milp.jl
@@ -211,7 +211,7 @@ function constraint_pump_head_gain(wm::AbstractMILPModel, n::Int, a::Int, node_f
     c_4 = JuMP.@constraint(wm.model, lambda[a, end] <= x_pw[a, end])
 
     # Generate a set of uniform flow breakpoints.
-    breakpoints = range(_q_eps, stop=JuMP.upper_bound(q), length=pump_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)
@@ -357,7 +357,7 @@ function objective_owf(wm::AbstractMILPModel)
                 q_ub = JuMP.upper_bound(q)
 
                 # Generate a set of uniform flow and cubic function breakpoints.
-                breakpoints = range(_q_eps, stop=q_ub, length=pump_breakpoints)
+                breakpoints = range(0.0, stop=q_ub, length=pump_breakpoints)
                 f = _calc_cubic_flow_values(collect(breakpoints), curve_fun)
 
                 # Get pump efficiency data.
diff --git a/src/form/milpr.jl b/src/form/milpr.jl
index 845b0fa5..254737b8 100644
--- a/src/form/milpr.jl
+++ b/src/form/milpr.jl
@@ -13,7 +13,7 @@ function _get_head_loss_oa(q::JuMP.VariableRef, q_hat::Float64, alpha::Float64)
 end
 
 
-function _get_head_loss_cv_oa(q::JuMP.VariableRef, z::Union{JuMP.VariableRef, JuMP.GenericAffExpr}, q_hat::Float64, alpha::Float64)
+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)
 end
 
@@ -42,7 +42,7 @@ function constraint_pump_head_gain(wm::AbstractMILPRModel, n::Int, a::Int, node_
     c_4 = JuMP.@constraint(wm.model, lambda[a, end] <= x_pw[a, end])
 
     # Add a constraint for the flow piecewise approximation.
-    breakpoints = range(_q_eps, stop=JuMP.upper_bound(qp), length=pump_breakpoints)
+    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)
 
@@ -99,12 +99,15 @@ function constraint_pipe_head_loss(wm::AbstractMILPRModel, n::Int, a::Int, node_
     pipe_breakpoints = get(wm.ext, :pipe_breakpoints, 0)
     if pipe_breakpoints <= 0 return end
 
+    # 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)
 
     for qp_hat in range(0.0, stop=JuMP.upper_bound(qp), length=pipe_breakpoints)
-        lhs = r * _get_head_loss_oa(qp, qp_hat, alpha)
-        c_p = JuMP.@constraint(wm.model, lhs <= inv(L) * dhp)
+        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
 
@@ -112,8 +115,8 @@ function constraint_pipe_head_loss(wm::AbstractMILPRModel, n::Int, a::Int, node_
     qn, dhn = var(wm, n, :qn_pipe, a), var(wm, n, :dhn_pipe, a)
 
     for qn_hat in range(0.0, stop=JuMP.upper_bound(qn), length=pipe_breakpoints)
-        lhs = r * _get_head_loss_oa(qn, qn_hat, alpha)
-        c_n = JuMP.@constraint(wm.model, lhs <= inv(L) * dhn)
+        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
 end
@@ -125,14 +128,16 @@ function constraint_check_valve_head_loss(wm::AbstractMILPRModel, n::Int, a::Int
     if pipe_breakpoints <= 0 return end
 
     # Get common variables for outer approximation constraints.
-    qp, z = var(wm, n, :qp_check_valve, a), var(wm, n, :z_check_valve, a)
-    dhp = var(wm, n, :dhp_check_valve, a)
+    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 = _get_head_loss_cv_oa(qp, z, qp_hat, ref(wm, n, :alpha))
-        c_p = JuMP.@constraint(wm.model, r * lhs <= inv(L) * dhp)
-        append!(con(wm, n, :head_loss)[a], [c_p])
+        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
 
@@ -145,20 +150,24 @@ function constraint_shutoff_valve_head_loss(wm::AbstractMILPRModel, n::Int, a::I
     # 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 = var(wm, n, :y_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 = _get_head_loss_cv_oa(qp, y, qp_hat, ref(wm, n, :alpha))
-        c_p = JuMP.@constraint(wm.model, L * r * lhs <= dhp)
-        append!(con(wm, n, :head_loss)[a], [c_p])
+        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])
     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 = _get_head_loss_cv_oa(qn, (1.0 - y), qn_hat, ref(wm, n, :alpha))
-        c_n = JuMP.@constraint(wm.model, L * r * lhs <= dhn)
-        append!(con(wm, n, :head_loss)[a], [c_n])
+        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])
     end
 end
 
@@ -192,7 +201,7 @@ function objective_owf(wm::AbstractMILPRModel)
                 qp_ub = JuMP.upper_bound(qp)
 
                 # Generate a set of uniform flow and cubic function breakpoints.
-                breakpoints = range(_q_eps, stop=qp_ub, length=pump_breakpoints)
+                breakpoints = range(0.0, stop=qp_ub, length=pump_breakpoints)
                 f = _calc_cubic_flow_values(collect(breakpoints), curve_fun)
 
                 # Get pump efficiency data.
diff --git a/src/form/nlp.jl b/src/form/nlp.jl
index 7da9337a..27e0088f 100644
--- a/src/form/nlp.jl
+++ b/src/form/nlp.jl
@@ -23,7 +23,7 @@ function constraint_check_valve_head_loss(wm::AbstractNLPModel, n::Int, a::Int,
     # 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 = JuMP.lower_bound(h_i) - JuMP.upper_bound(h_j)
+    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
@@ -44,8 +44,8 @@ function constraint_shutoff_valve_head_loss(wm::AbstractNLPModel, n::Int, a::Int
     # 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 = JuMP.lower_bound(h_i) - JuMP.upper_bound(h_j)
-    dh_ub = JuMP.upper_bound(h_i) - JuMP.lower_bound(h_j)
+    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
diff --git a/src/form/undirected_flow.jl b/src/form/undirected_flow.jl
index 9ce92ac7..68d66cfc 100644
--- a/src/form/undirected_flow.jl
+++ b/src/form/undirected_flow.jl
@@ -97,23 +97,22 @@ function constraint_check_valve_common(wm::AbstractUndirectedModel, n::Int, a::I
     # Get flow and check valve status variables.
     q, z = var(wm, n, :q_check_valve, a), var(wm, n, :z_check_valve, a)
 
-    # If the check valve is open, flow must be appreciably nonnegative.
+    # If the check valve is closed, flow 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 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 check valve is open, negative head loss is not possible.
     dh_lb = JuMP.lower_bound(h_i) - JuMP.upper_bound(h_j)
-    c_3 = JuMP.@constraint(wm.model, h_i - h_j >= (1.0 - z) * dh_lb)
+    c_2 = JuMP.@constraint(wm.model, h_i - h_j >= (1.0 - z) * min(0.0, dh_lb))
 
     # When the check valve is closed, positive head loss is not possible.
     dh_ub = JuMP.upper_bound(h_i) - JuMP.lower_bound(h_j)
-    c_4 = JuMP.@constraint(wm.model, h_i - h_j <= z * dh_ub)
+    c_3 = JuMP.@constraint(wm.model, h_i - h_j <= z * dh_ub)
 
     # Append the constraint array.
-    append!(con(wm, n, :check_valve, a), [c_1, c_2, c_3, c_4])
+    append!(con(wm, n, :check_valve, a), [c_1, c_2, c_3])
 end
 
 
@@ -168,9 +167,9 @@ function constraint_pump_common(wm::AbstractUndirectedModel, n::Int, a::Int, nod
     c_2 = JuMP.@constraint(wm.model, q >= _q_eps * z)
 
     # If the pump is off, decouple the head difference relationship.
-    dhn_lb = JuMP.lower_bound(h_j) - JuMP.upper_bound(h_i)
+    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 = JuMP.upper_bound(h_j) - JuMP.lower_bound(h_i)
+    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))
 
     # Append the constraint array.
@@ -181,6 +180,15 @@ function constraint_pipe_common(wm::AbstractUndirectedModel, n::Int, a::Int, nod
     # For undirected formulations, there are no constraints, here.
 end
 
+function constraint_node_directionality(
+    wm::AbstractWaterModel, n::Int, i::Int, check_valve_fr::Array{Int},
+    check_valve_to::Array{Int}, pipe_fr::Array{Int}, pipe_to::Array{Int},
+    pump_fr::Array{Int}, pump_to::Array{Int}, pressure_reducing_valve_fr::Array{Int},
+    pressure_reducing_valve_to::Array{Int}, shutoff_valve_fr::Array{Int},
+    shutoff_valve_to::Array{Int})
+    # For undirected formulations, there are no constraints, here.
+end
+
 function constraint_sink_directionality(
     wm::AbstractWaterModel, n::Int, i::Int, check_valve_fr::Array{Int},
     check_valve_to::Array{Int}, pipe_fr::Array{Int}, pipe_to::Array{Int},
diff --git a/src/io/epanet.jl b/src/io/epanet.jl
index de6fc83d..98a868f1 100644
--- a/src/io/epanet.jl
+++ b/src/io/epanet.jl
@@ -285,7 +285,7 @@ function _update_reservoir_ts!(data::Dict{String,<:Any})
     for (i, reservoir) in data["time_series"]["reservoir"]
         key = findfirst(x -> x["source_id"][2] == i, data["reservoir"])
         node_index = data["reservoir"][key]["node"]
-        node_ts[key] = reservoir
+        node_ts[string(node_index)] = reservoir
     end
 
     # Update the reservoir entry in the time series dictionary.
@@ -671,6 +671,7 @@ function _read_junction!(data::Dict{String,<:Any})
         # 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
 
         # Parse the elevation of the junction (in meters).
         if data["flow_units"] == "LPS" || data["flow_units"] == "CMH"
@@ -994,6 +995,7 @@ function _read_reservoir!(data::Dict{String,<:Any})
         # Initialize the reservoir entry to be added.
         reservoir = Dict{String,Any}("name" => current[1], "status" => 1)
         reservoir["source_id"] = ["reservoir", current[1]]
+        reservoir["dispatchable"] = false
 
         # Parse the head of the reservoir (in meters).
         if data["flow_units"] == "LPS" || data["flow_units"] == "CMH"
@@ -1049,6 +1051,7 @@ function _read_tank!(data::Dict{String,<:Any})
         # Initialize the tank entry to be added.
         tank = Dict{String,Any}("name" => current[1], "status" => 1)
         tank["source_id"] = ["tank", current[1]]
+        tank["dispatchable"] = false
 
         # Store all measurements associated with tanks in metric units.
         if data["flow_units"] == "LPS" || data["flow_units"] == "CMH"
@@ -1126,11 +1129,13 @@ end
 
 
 function _read_valve!(data::Dict{String,<:Any})
-    # Create a shorthand variable for the pressure reducing valve type.
+    # Create a shorthand variable for each of the valve types.
     prv_type = "pressure_reducing_valve"
+    tcv_type = "throttle_control_valve"
 
     # Initialize dictionaries associated with valves.
     data[prv_type] = Dict{String,Dict{String,Any}}()
+    data[tcv_type] = Dict{String,Dict{String,Any}}()
 
     # Initialize a temporary index to be updated while parsing.
     index::Int64 = 0
@@ -1148,6 +1153,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]]
+        elseif uppercase(current[5]) == "TCV"
+            valve["source_id"] = [tcv_type, 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.")
         end
diff --git a/src/prob/des.jl b/src/prob/des.jl
index bc48ca8d..676296d9 100644
--- a/src/prob/des.jl
+++ b/src/prob/des.jl
@@ -38,6 +38,7 @@ function build_des(wm::AbstractWaterModel)
     # Flow conservation at all nodes.
     for (i, node) in ref(wm, :node)
         constraint_flow_conservation(wm, i)
+        constraint_node_directionality(wm, i)
     end
 
     # Set source node hydraulic heads.
diff --git a/src/prob/owf.jl b/src/prob/owf.jl
index 55c04199..01690348 100644
--- a/src/prob/owf.jl
+++ b/src/prob/owf.jl
@@ -22,7 +22,7 @@ function build_mn_owf(wm::AbstractWaterModel)
     network_ids = sort(collect(nw_ids(wm)))
 
     # Ensure tanks recover their initial volume.
-    n_1, n_f = [network_ids[1], network_ids[end]]
+    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)
diff --git a/src/prob/wf.jl b/src/prob/wf.jl
index af71cbc9..eea204c8 100644
--- a/src/prob/wf.jl
+++ b/src/prob/wf.jl
@@ -18,6 +18,7 @@ function build_wf(wm::AbstractWaterModel)
     variable_shutoff_valve_indicator(wm)
 
     # Component-specific variables.
+    variable_dispatchable_junction(wm)
     variable_pump_operation(wm)
     variable_reservoir(wm)
     variable_tank(wm)
@@ -25,6 +26,7 @@ function build_wf(wm::AbstractWaterModel)
     # Flow conservation at all nodes.
     for (i, node) in ref(wm, :node)
         constraint_flow_conservation(wm, i)
+        constraint_node_directionality(wm, i)
     end
 
     # Head loss along pipes without valves.
@@ -60,9 +62,13 @@ function build_wf(wm::AbstractWaterModel)
 
     # Constrain flow directions based on demand.
     for (i, junction) in ref(wm, :junction)
-        if junction["demand"] > 0.0
+        if !junction["dispatchable"] && junction["demand"] > 0.0
             constraint_sink_directionality(wm, junction["node"])
-        elseif junction["demand"] < 0.0
+        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
     end
@@ -97,6 +103,7 @@ function build_mn_wf(wm::AbstractWaterModel)
         variable_pump_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)
@@ -104,6 +111,7 @@ function build_mn_wf(wm::AbstractWaterModel)
         # Flow conservation at all nodes.
         for i in ids(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.
@@ -139,9 +147,9 @@ function build_mn_wf(wm::AbstractWaterModel)
 
         # Constrain flow directions based on demand.
         for (i, junction) in ref(wm, :junction; nw=n)
-            if junction["demand"] > 0.0
+            if !junction["dispatchable"] && junction["demand"] > 0.0
                 constraint_sink_directionality(wm, junction["node"]; nw=n)
-            elseif junction["demand"] < 0.0
+            elseif !junction["dispatchable"] && junction["demand"] < 0.0
                 constraint_source_directionality(wm, junction["node"]; nw=n)
             end
         end
diff --git a/src/util/obbt.jl b/src/util/obbt.jl
new file mode 100644
index 00000000..b9033047
--- /dev/null
+++ b/src/util/obbt.jl
@@ -0,0 +1,400 @@
+### Optimization-based Bound Tightening for Optimal Water Flow Problems ###
+### This is built upon https://github.com/lanl-ansi/PowerModels.jl/blob/master/src/util/obbt.jl
+
+"""
+Iteratively tighten bounds on head and flow variables.
+
+The function can be invoked on any convex relaxation which explicitly has these variables.
+By default, the function uses the MICP-R relaxation for performing bound tightening.
+
+# Example
+
+The function can be invoked as follows:
+
+```
+ipopt = JuMP.optimizer_with_attributes(Ipopt.Optimizer)
+run_obbt_owf!("examples/data/epanet/van_zyl.inp", ipopt)
+```
+
+# Keyword Arguments
+* `model_type`: relaxation to use for performing bound tightening.
+* `max_iter`: maximum number of bound tightening iterations to perform.
+* `min_width`: domain width beyond which bound tightening is not performed.
+* `precision`: decimal precision to round the tightened bounds to.
+* `time_limit`: maximum amount of time (seconds) for the bound tightening algorithm.
+* `upper_bound`: can be used to specify a local feasible solution objective for the owf problem.
+* `upper_bound_constraint`: boolean option that can be used to add an additional
+   constraint to reduce the search space of each of the bound tightening
+   solves. This cannot be set to `true` without specifying an upper bound.
+* `use_reduced_network`: boolean option that specifies whether or not to use a reduced,
+   snapshot, dispatchable version of the origin multinetwork problem for bound tightening.
+"""
+function run_obbt_owf!(file::String, optimizer; kwargs...)
+    data = WaterModels.parse_file(file)
+    return run_obbt_owf!(data, optimizer; kwargs...)
+end
+
+
+function _relax_model!(wm::AbstractWaterModel)
+    # Further relax all binary variables in the model.
+    bin_vars = filter(v -> JuMP.is_binary(v), JuMP.all_variables(wm.model))
+    JuMP.unset_binary.(bin_vars) # Make all binary variables free and continuous.
+    JuMP.set_lower_bound.(bin_vars, 0.0) # Lower-bound the relaxed binary variables.
+    JuMP.set_upper_bound.(bin_vars, 1.0) # Upper-bound the relaxed binary variables.
+end
+
+
+function _get_node_bound_dict(wm::AbstractWaterModel, nw::Int)
+    return Dict(string(i) => Dict("h_min"=>JuMP.lower_bound(var(wm, nw, :h, i)),
+        "h_max"=>JuMP.upper_bound(var(wm, nw, :h, i))) for i in ids(wm, nw, :node))
+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))
+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))
+        data[string(type)] = _get_edge_bound_dict(wm, wm.cnw, type)
+    end
+
+    return data
+end
+
+
+function _create_modifications_mn(wm::AbstractWaterModel)
+    data = Dict{String,Any}("nw"=>Dict{String,Any}(), "per_unit"=>false, "multinetwork"=>true)
+
+    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))
+            dnw[string(type)] = _get_edge_bound_dict(wm, nw, type)
+        end
+    end
+
+    return data
+end
+
+
+function _get_head_index_set(wm::AbstractWaterModel; width::Float64=1.0e-3, prec::Float64=1.0e-4)
+    return vcat([[(nw, :node, :h, i, width, prec) for i in ids(wm, nw, :node)]
+        for nw in sort(collect(nw_ids(wm)))]...)
+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]
+    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_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])
+
+    if "nw" in keys(data)
+        data_index = data["nw"][nw_str][comp_str][index_str]
+    else
+        data_index = data[comp_str][index_str]
+    end
+
+    return data_index[var_str * "_min"], data_index[var_str * "_max"]
+end
+
+
+function _update_modifications!(data::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
+
+        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]
+        else
+            data[comp_c][index_str][var_str * "_min"] = bounds[i][1]
+            data[comp_c][index_str][var_str * "_max"] = bounds[i][2]
+        end
+    end
+end
+
+
+function _get_average_widths(data::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"]
+        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))
+        end
+    end
+
+    Memento.info(_LOGGER, message[1:end-2] * ".")
+    return avg_vals
+end
+
+
+function _get_average_widths_mn(data::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"]
+        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"])
+            message *= "q_$(type) -> $(q * inv(q_length)), "
+            avg_vals = vcat(avg_vals, q * inv(q_length))
+        end
+    end
+
+    Memento.info(_LOGGER, message[1:end-2] * ".")
+    return avg_vals
+end
+
+
+"Translate a multinetwork dataset to a snapshot dataset with dispatchable components."
+function _make_reduced_data!(ts_data::Dict{String,<:Any})
+    for comp_type in keys(ts_data["time_series"])
+        # If not a component type (Dict), skip parsing.
+        !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"])
+            elseif comp_type == "node"
+                ts_data["node"][i]["h_min"] = minimum(comp["head"])
+                ts_data["node"][i]["h_max"] = maximum(comp["head"])
+            end
+        end
+    end
+
+    for (i, tank) in ts_data["tank"]
+        tank["dispatchable"] = true
+    end
+
+    for (i, reservoir) in ts_data["reservoir"]
+        if "node" in keys(ts_data["time_series"])
+            if string(reservoir["node"]) in keys(ts_data["time_series"]["node"])
+                reservoir["dispatchable"] = true
+            end
+        end
+    end
+end
+
+
+function _revert_reduced_data!(ts_data::Dict{String,<:Any})
+    for comp_type in keys(ts_data["time_series"])
+        # If not a component type (Dict), skip parsing.
+        !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")
+            end
+        end
+    end
+
+    for (i, tank) in ts_data["tank"]
+        tank["dispatchable"] = false
+    end
+
+    for (i, reservoir) in ts_data["reservoir"]
+        node_id = string(reservoir["node"])
+        reservoir["dispatchable"] = false
+        delete!(ts_data["node"][node_id], "h_min")
+        delete!(ts_data["node"][node_id], "h_max")
+    end
+end
+
+
+function run_obbt_owf!(data::Dict{String,<:Any}, optimizer; use_reduced_network::Bool = true,
+    model_type::Type = MICPRWaterModel, 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...)
+    # 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.")
+
+    if !_IM.ismultinetwork(data) && use_reduced_network
+        _make_reduced_data!(data)
+        build_type = WaterModels.build_wf
+    elseif _IM.ismultinetwork(data) && !use_reduced_network
+        build_type = WaterModels.build_mn_owf
+    else
+        build_type = WaterModels.build_mn_owf
+        Memento.error(_LOGGER, "[OBBT] Ensure input network is the correct type.")
+    end
+
+    # Check for keyword argument inconsistencies.
+    _check_obbt_options(upper_bound, upper_bound_constraint)
+
+    # Instantiate the bound tightening model and relax integrality, if specified.
+    bt = instantiate_model(data, model_type, build_type; ext=ext)
+    upper_bound_constraint && _constraint_obj_bound(bt, upper_bound)
+    relaxed && _relax_model!(bt) # Relax integrality, if required.
+
+    # Build the dictionary and sets that will store to the network.
+    modifications = use_reduced_network ? _create_modifications_reduced(bt) : _create_modifications_mn(bt)
+    head_index_set = _get_head_index_set(bt; width=min_width, prec=precision)
+    flow_index_set = _get_flow_index_set(bt; width=min_width, prec=precision)
+    var_index_set = vcat(head_index_set, flow_index_set)
+    bounds = [_get_existing_bounds(modifications, vid) for vid in var_index_set]
+
+    # Get the initial average bound widths.
+    if use_reduced_network
+        avg_widths_initial = _get_average_widths(modifications)
+    else
+        avg_widths_initial = _get_average_widths_mn(modifications)
+    end
+
+    # Initialize algorithm termination metadata.
+    current_iteration, time_elapsed = 1, 0.0
+    parallel_time_elapsed, terminate = 0.0, false
+
+    # Set the optimizer for the bound tightening model.
+    JuMP.set_optimizer(bt.model, optimizer)
+
+    while !terminate # Algorithmic loop.
+        Memento.info(_LOGGER, "[OBBT] Starting iteration $(current_iteration).")
+
+        # Loop over variables and tighten bounds.
+        for (i, id) in enumerate(var_index_set)
+            # If current bounds are within the minimum variable width, skip.
+            bounds[i][2] - bounds[i][1] < id[5] && continue
+
+            # Perform optimization-based bound tightening for lower and upper bounds.
+            time_elapsed += @elapsed lb = _solve_bound(bt, id, _MOI.MIN_SENSE, bounds[i][1])
+            time_elapsed > time_limit && ((terminate = true) && break)
+            time_elapsed += @elapsed ub = _solve_bound(bt, id, _MOI.MAX_SENSE, bounds[i][2])
+            time_elapsed > time_limit && ((terminate = true) && break)
+
+            # Compute corrected bounds from the optimization-based bounds.
+            bounds[i] = _compute_corrected_bounds(lb, ub, bounds[i][1], bounds[i][2], id[5], id[6])
+        end
+
+        # Update the iteration counter.
+        current_iteration += 1
+
+        # Set the termination variable if max iterations is exceeded.
+        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)
+        _IM.update_data!(data, modifications)
+
+        # Get the current average bound widths.
+        if use_reduced_network
+            avg_widths_final = _get_average_widths(modifications)
+        else
+            avg_widths_final = _get_average_widths_mn(modifications)
+        end
+
+        if maximum(avg_widths_initial - avg_widths_final) < improvement_tol
+            terminate = true
+        end
+
+        if !terminate
+            # Set the initial average widths to the last average widths.
+            avg_widths_initial = avg_widths_final
+
+            # Instantiate the new model and set the optimizer.
+            bt = instantiate_model(data, model_type, build_type; ext=ext)
+            upper_bound_constraint && _constraint_obj_bound(bt, upper_bound)
+            relaxed && _relax_model!(bt)
+            JuMP.set_optimizer(bt.model, optimizer)
+        end
+    end
+
+    if use_reduced_network
+        _revert_reduced_data!(data)
+    end
+end
+
+
+function _check_obbt_options(ub::Float64, ub_constraint::Bool)
+    if ub_constraint && isinf(ub)
+        Memento.error(_LOGGER, "[OBBT] The option \"upper_bound_constraint\" cannot be set to true without specifying an upper bound.")
+    end
+end
+
+
+function _compute_corrected_bounds(lb::Float64, ub::Float64, lb_old::Float64, ub_old::Float64, width::Float64, prec::Float64)
+    # Convert new bounds to the desired decimal precision.
+    lb = max(floor(inv(prec) * lb) * prec, lb_old)
+    ub = min(ceil(inv(prec) * ub) * prec, ub_old)
+
+    # If the bounds satisfy the minimum width, return them.
+    ub - lb >= width && return lb, ub
+
+    # Compute the mean of the bounds.
+    mean = 0.5 * (lb + ub)
+
+    # Return bounds satisfying the minimum width.
+    if mean - 0.5 * width < lb_old
+        return lb_old, ub_old + width
+    elseif mean + 0.5 * width > ub_old
+        return ub_old - width, ub_old
+    else
+        return mean - 0.5 * width, mean + 0.5 * width
+    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])
+
+    # Optimize the variable (or affine expression) being tightened.
+    JuMP.@objective(wm.model, sense, variable)
+    JuMP.optimize!(wm.model)
+
+    # 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.
+        candidate = JuMP.objective_value(wm.model)
+        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)
+        return start # Optimization was not successful. Return the starting bound.
+    end
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index 9368875b..2492efd4 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -44,6 +44,8 @@ include("common.jl")
 
     include("des.jl")
 
+    include("util.jl")
+
     include("warm_start.jl")
 
 end
diff --git a/test/util.jl b/test/util.jl
new file mode 100644
index 00000000..38008e00
--- /dev/null
+++ b/test/util.jl
@@ -0,0 +1,14 @@
+@testset "Optimization-based Bound Tightening (Reduced Network)" begin
+    network = WaterModels.parse_file("../test/data/epanet/multinetwork/pump-hw-lps.inp")
+    run_obbt_owf!(network, ipopt)
+    @test haskey(network["pump"]["1"], "q_min")
+    @test haskey(network["pump"]["1"], "q_max")
+end
+
+@testset "Optimization-based Bound Tightening (Multinetwork)" begin
+    network = WaterModels.parse_file("../test/data/epanet/multinetwork/pump-hw-lps.inp")
+    network_mn = WaterModels.make_multinetwork(network)
+    run_obbt_owf!(network_mn, ipopt, use_reduced_network=false)
+    @test haskey(network_mn["nw"]["1"]["pump"]["1"], "q_min")
+    @test haskey(network_mn["nw"]["1"]["pump"]["1"], "q_max")
+end
diff --git a/test/wf.jl b/test/wf.jl
index 4093d1d8..6a334aac 100644
--- a/test/wf.jl
+++ b/test/wf.jl
@@ -306,7 +306,6 @@ end
         @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
-        @test result["solution"]["nw"]["3"]["pump"]["1"]["g"] <= 99.60
 
         result = run_mn_wf(network, MILPWaterModel, cbc, ext=Dict(:pump_breakpoints=>5))
         @test result["termination_status"] == OPTIMAL