From bc82dd40f5ba04299f7509857a9804244e8c8005 Mon Sep 17 00:00:00 2001 From: Byron Tasseff Date: Thu, 30 Jul 2020 15:26:37 -0600 Subject: [PATCH 01/24] Tidying up constraints involving head difference bounds. --- src/form/micp.jl | 5 ++--- src/form/milpr.jl | 22 +++++++++++++--------- src/form/nlp.jl | 6 +++--- src/form/undirected_flow.jl | 6 +++--- 4 files changed, 21 insertions(+), 18 deletions(-) diff --git a/src/form/micp.jl b/src/form/micp.jl index 5b1b560a..44256967 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]) @@ -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) diff --git a/src/form/milpr.jl b/src/form/milpr.jl index 845b0fa5..963ff302 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 @@ -130,7 +130,7 @@ function constraint_check_valve_head_loss(wm::AbstractMILPRModel, n::Int, a::Int # 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)) + lhs = _get_head_loss_oa_z(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]) end @@ -145,20 +145,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 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..a74bfe39 100644 --- a/src/form/undirected_flow.jl +++ b/src/form/undirected_flow.jl @@ -106,7 +106,7 @@ function constraint_check_valve_common(wm::AbstractUndirectedModel, n::Int, a::I # 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_3 = 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) @@ -168,9 +168,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. From 82b0548a6e44279c6941581f66a73d59347ab88f Mon Sep 17 00:00:00 2001 From: Byron Tasseff Date: Mon, 3 Aug 2020 10:23:30 -0600 Subject: [PATCH 02/24] Making progress on OBBT utility. --- .travis.yml | 4 +- src/WaterModels.jl | 1 + src/util/obbt.jl | 369 +++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 372 insertions(+), 2 deletions(-) create mode 100644 src/util/obbt.jl diff --git a/.travis.yml b/.travis.yml index 41f1cd5c..ed6b4084 100644 --- a/.travis.yml +++ b/.travis.yml @@ -4,7 +4,7 @@ os: - osx julia: - 1.0 - - 1.4 + - 1.5 - nightly codecov: true jobs: @@ -12,7 +12,7 @@ jobs: - julia: nightly include: - stage: "Documentation" - julia: 1.4 + julia: 1.5 os: linux script: - julia --project=docs/ -e 'using Pkg; Pkg.instantiate(); Pkg.develop(PackageSpec(path=pwd()))' 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/util/obbt.jl b/src/util/obbt.jl new file mode 100644 index 00000000..40f07d28 --- /dev/null +++ b/src/util/obbt.jl @@ -0,0 +1,369 @@ +### Optimization-based Bound Tightening for Optimal Water Flow Problems ### +### This is built upon https://github.com/lanl-ansi/WaterModels.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. +Interested readers are referred to the paper "Strengthening Convex Relaxations with Bound +Tightening for Water Network Optimization." + +# Example + +The function can be invoked as follows: + +``` +ipopt = JuMP.optimizer_with_attributes(Ipopt.Optimizer) +data, stats = run_obbt_owf!("examples/data/epanet/van_zyl.inp", ipopt) +``` + +`data` contains the parsed network data with tightened bounds. +`stats` contains information output from the bound tightening algorithm. + +# Keyword Arguments +* `model_type`: relaxation to use for performing bound tightening. +* `max_iter`: maximum number of bound tightening iterations to perform. +* `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. +* `rel_gap_tol`: tolerance used to terminate the algorithm when the objective + value of the relaxation is close to the upper bound specified using the + `upper_bound` keyword. +* `min_bound_width`: domain beyond which bound tightening is not performed. +* `termination`: Bound-tightening algorithm terminates if the improvement in + the average or maximum bound improvement, specified using either the + `termination = :avg` or the `termination = :max` option, is less than + `improvement_tol`. +* `precision`: number of decimal digits to round the tightened bounds to. +""" +function run_obbt_owf!(file::String, optimizer; kwargs...) + data = WaterModels.parse_file(file) + return run_obbt_owf!(data, optimizer; kwargs...) +end + +function _get_obbt_node_var_lb(wm::AbstractWaterModel, v::Symbol) + network_ids = sort(collect(nw_ids(wm))) + return Dict(nw => Dict(i => JuMP.lower_bound(var(wm, nw, v, i)) + for i in ids(wm, nw, :node)) for nw in network_ids) +end + +function _get_obbt_node_var_ub(wm::AbstractWaterModel, v::Symbol) + network_ids = sort(collect(nw_ids(wm))) + return Dict(nw => Dict(i => JuMP.upper_bound(var(wm, nw, v, i)) + for i in ids(wm, nw, :node)) for nw in network_ids) +end + +function run_obbt_owf!(data::Dict{String,<:Any}, optimizer; + model_type::Type = MICPRWaterModel, + max_iter::Int = 100, + time_limit::Float64 = 3600.0, + upper_bound::Float64 = Inf, + upper_bound_constraint::Bool = false, + rel_gap_tol::Float64 = Inf, + min_bound_width::Float64 = 1.0e-2, + improvement_tol::Float64 = 1.0e-3, + precision::Int = 4, + termination::Symbol = :avg, + kwargs...) + + Memento.info(_LOGGER, "Maximum number of OBBT iterations set to default value of $(max_iter).") + Memento.info(_LOGGER, "Maximum time limit for OBBT set to default value of $(time_limit) seconds.") + model_relaxation = instantiate_model(data, model_type, WaterModels.build_mn_owf) + + # Check for model_type compatability with OBBT. + _check_variables(model_relaxation) + + # Check for other keyword argument consistencies. + _check_obbt_options(upper_bound, rel_gap_tol, upper_bound_constraint) + + # Check termination norm criteria for OBBT. + if termination != :avg && termination != :max + Memento.error(_LOGGER, "OBBT termination criteria can only be :max or :avg.") + end + + # Further relax all binary variables in the model. + bin_vars = filter(v -> JuMP.is_binary(v), JuMP.all_variables(model_relaxation.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. + + # Declare possible pass statuses. + pass_statuses = [_MOI.LOCALLY_SOLVED, _MOI.OPTIMAL] + + # Compute initial relative gap between relaxation objective and upper_bound. + result_relaxation = optimize_model!(model_relaxation, optimizer=optimizer) + current_relaxation_objective = result_relaxation["objective"] + + if upper_bound < current_relaxation_objective + Memento.error(_LOGGER, "The upper bound provided to OBBT is not a valid OWF upper bound.") + end + + if !(result_relaxation["termination_status"] in pass_statuses) + termination_status = result_relaxation["termination_status"] + Memento.warn(_LOGGER, "Initial relaxation solve status is $(termination_status).") + + if termination_status == :SubOptimal + Memento.warn(_LOGGER, "Continuing with the bound tightening algorithm.") + end + end + + current_rel_gap = Inf + + if !isinf(upper_bound) + current_rel_gap = (upper_bound - current_relaxation_objective) * inv(upper_bound) + Memento.info(_LOGGER, "Initial relaxation gap is $(current_rel_gap).") + end + + model_bt = instantiate_model(data, model_type, WaterModels.build_mn_owf) + upper_bound_constraint && _constraint_obj_bound(model_bt, upper_bound) + + # Further relax all binary variables in the bound tightening model. + bin_vars = filter(v -> JuMP.is_binary(v), JuMP.all_variables(model_bt.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. + + # Initialize the statistics dictionary. + stats = Dict{String,Any}() + stats["model_type"] = model_type + stats["initial_relaxation_objective"] = current_relaxation_objective + stats["initial_rel_gap_from_ub"] = current_rel_gap + stats["upper_bound"] = upper_bound + + h_lb = _get_obbt_node_var_lb(model_bt, :h) + h_ub = _get_obbt_node_var_ub(model_bt, :h) + stats["h_range_init"], stats["avg_h_range_init"] = 0.0, 0.0 + + for nw in sort(collect(nw_ids(model_bt))) + node_ids = ids(model_bt, nw, :node) + h_range_nw = sum(h_ub[nw][i] - h_lb[nw][i] for i in node_ids) + stats["h_range_init"] += h_range_nw + stats["avg_h_range_init"] += h_range_nw * inv(length(node_ids)) + end + + # Initialize algorithm termination metadata. + h_range_final, total_h_reduction = 0.0, Inf + max_h_reduction, avg_h_reduction = Inf, Inf + final_relaxation_objective = NaN + current_iteration, time_elapsed = 0, 0.0 + parallel_time_elapsed, terminate = 0.0, false + + # Set whether or not the algorithm should immediately terminate. + if termination == :avg + terminate = avg_h_reduction <= improvement_tol && avg_q_reduction <= improvement_tol + elseif termination == :max + terminate = max_h_reduction <= improvement_tol && max_q_reduction <= improvement_tol + end + + while !terminate # Algorithmic loop. + # Set important metadata describing the current round. + iter_start_time, max_h_iteration_time = time(), 0.0 + total_h_reduction, avg_h_reduction, max_h_reduction = 0.0, 0.0, 0.0 + + # Loop over all subnetworks in the multinetwork. + for nw in sort(collect(nw_ids(model_bt))) + + # Loop over all nodes in the network. + for i in ids(model_bt, nw, :node) + println(nw, " ", i) + + # If the lower and upper bounds are already close, skip. + if h_ub[nw][i] - h_lb[nw][i] < min_bound_width + continue + end + + # Store metadata for variable tightening. + lb, ub, start_time = NaN, NaN, time() + + # Minimize the variable whose bounds are being tightened. + JuMP.@objective(model_bt.model, _MOI.MIN_SENSE, var(model_bt, nw, :h, i)) + result_bt = optimize_model!(model_bt, optimizer=optimizer) + + # Store the new lower bound or print a warning and continue. + if result_bt["termination_status"] in pass_statuses + # Compute the candidate lower bound at the specified precision. + val = JuMP.objective_value(model_bt.model) + lb_new = floor(10.0^precision * val) * inv(10.0^precision) + lb_new > h_lb[nw][i] && (lb = lb_new) # Store the new lower bound. + else + Memento.warn(_LOGGER, "BT minimization for h[$(i)] errored. Adjust tolerances.") + continue + end + + # Maximize the variable whose bounds are being tightened. + JuMP.@objective(model_bt.model, _MOI.MAX_SENSE, var(model_bt, nw, :h, i)) + result_bt = optimize_model!(model_bt, optimizer=optimizer) + + # Store the new lower bound or print a warning and continue. + if result_bt["termination_status"] in pass_statuses + # Compute the candidate lower bound at the specified precision. + val = JuMP.objective_value(model_bt.model) + ub_new = ceil(10.0^precision * val) * inv(10.0^precision) + ub_new < h_ub[nw][i] && (ub = ub_new) # Store the new lower bound. + else + Memento.warn(_LOGGER, "BT maximization for h[$(i)] errored. Adjust tolerances.") + continue + end + + # Compute the time required for the variable's bound tightening. + end_time = time() - start_time + max_h_iteration_time = max(end_time, max_h_iteration_time) + + # Perform sanity checks on the new bounds. + if lb > ub # Check if the lower bound exceeds the upper bound. + Memento.warn(_LOGGER, "BT lb > ub for h[$(i)]. Adjust tolerances.") + continue + end + + # Revert to old bounds if new bounds are nonsensical or NaN. + (!isnan(lb) && lb > h_ub[i]) && (lb = h_lb[i]) + (!isnan(ub) && ub < h_lb[i]) && (ub = h_ub[i]) + isnan(lb) && (lb = h_lb[i]) + isnan(ub) && (ub = h_ub[i]) + end + end + end + + # # vm bound-reduction computation + # h_reduction = 0.0 + # if (ub - lb >= min_bound_width) + # h_reduction = (h_ub[node] - h_lb[node]) - (ub - lb) + # h_lb[node] = lb + # h_ub[node] = ub + # else + # mean = 0.5 * (ub + lb) + # if (mean - 0.5 * min_bound_width < h_lb[node]) + # lb = h_lb[node] + # ub = h_lb[node] + min_bound_width + # elseif (mean + 0.5 * min_bound_width > h_ub[node]) + # ub = h_ub[node] + # lb = h_ub[node] - min_bound_width + # else + # lb = mean - 0.5 * min_bound_width + # ub = mean + 0.5 * min_bound_width + # end + # h_reduction = (h_ub[node] - h_lb[node]) - (ub - lb) + # h_lb[node] = lb + # h_ub[node] = ub + # end + + # total_h_reduction += (h_reduction) + # max_h_reduction = max(h_reduction, max_h_reduction) + # end + # avg_h_reduction = total_h_reduction/length(nodes) + + # h_range_final = sum([h_ub[node] - h_lb[node] for node in nodes]) + + # q_range_final = sum([q_ub[bp] - q_lb[bp] for bp in pipes]) + + # parallel_time_elapsed += max(max_h_iteration_time, max_q_iteration_time) + + # time_elapsed += (time() - iter_start_time) + + # # populate the modifications, update the data, and rebuild the bound tightening model + # modifications = _create_modifications(model_bt, h_lb, h_ub, q_lb, q_ub) + # WaterModels.update_data!(data, modifications) + # model_bt = instantiate_model(data, model_type, WaterModels.build_mn_owf) + # (upper_bound_constraint) && (_constraint_obj_bound(model_bt, upper_bound)) + # vm = var(model_bt, :vm) + # q = var(model_bt, :q) + + # # run the qc relaxation for the updated bounds + # result_relaxation = run_owf(data, model_type::Type, optimizer) + + # if result_relaxation["termination_status"] in pass_statuses + # current_rel_gap = (upper_bound - result_relaxation["objective"])/upper_bound + # final_relaxation_objective = result_relaxation["objective"] + # else + # Memento.warn(_LOGGER, "relaxation solve failed in iteration $(current_iteration+1)") + # Memento.warn(_LOGGER, "using the previous iteration's gap to check relative gap stopping criteria") + # end + + # Memento.info(_LOGGER, "iteration $(current_iteration+1), vm range: $h_range_final, q range: $q_range_final, relaxation obj: $final_relaxation_objective") + + # # Update the termination flag. + # if termination == :avg + # terminate = avg_h_reduction <= improvement_tol && avg_q_reduction <= improvement_tol + # else + # terminate = max_h_reduction <= improvement_tol && max_q_reduction <= improvement_tol + # end + + # # iteration counter update + # current_iteration += 1 + # # check all the stopping criteria + # (current_iteration >= max_iter) && (Memento.info(_LOGGER, "maximum iteration limit reached"); break) + # (time_elapsed > time_limit) && (Memento.info(_LOGGER, "maximum time limit reached"); break) + # if (!isinf(rel_gap_tol)) && (current_rel_gap < rel_gap_tol) + # Memento.info(_LOGGER, "relative optimality gap < $rel_gap_tol") + # break + # end + + #end + + #pipes_vad_same_sign_count = 0 + #for (key, pipe) in data["pipe"] + # is_same_sign = (pipe["q_max"] >=0 && pipe["q_min"] >= 0) || (pipe["q_max"] <=0 && pipe["q_min"] <= 0) + # (is_same_sign) && (pipes_vad_same_sign_count += 1) + #end + + #stats["final_relaxation_objective"] = final_relaxation_objective + #stats["final_rel_gap_from_ub"] = isnan(upper_bound) ? Inf : current_rel_gap + #stats["h_range_final"] = h_range_final + #stats["avg_h_range_final"] = h_range_final * inv(length(nodes)) + + #stats["q_range_final"] = q_range_final + #stats["avg_q_range_final"] = q_range_final * inv(length(pipes)) + + #stats["run_time"] = time_elapsed + #stats["iteration_count"] = current_iteration + #stats["sim_parallel_run_time"] = parallel_time_elapsed + + #stats["vad_sign_determined"] = pipes_vad_same_sign_count + + return data, stats +end + + +function _check_variables(wm::AbstractWaterModel) + try + h = var(wm, :h) + catch err + (isa(error, KeyError)) && (Memento.error(_LOGGER, "OBBT is not supported for models without head variables.")) + end +end + + +function _check_obbt_options(ub::Float64, rel_gap::Float64, ub_constraint::Bool) + if ub_constraint && isinf(ub) + Memento.error(_LOGGER, "The option \"upper_bound_constraint\" cannot be set to true without specifying an upper bound.") + end + + if !isinf(rel_gap) && isinf(ub) + Memento.error(_LOGGER, "The option \"rel_gap_tol\" is specified without providing an upper bound.") + end +end + + +function _constraint_obj_bound(wm::AbstractWaterModel, bound) +end + + +function _create_modifications(wm::AbstractWaterModel, + h_lb::Dict{Any,Float64}, h_ub::Dict{Any,Float64}, + q_lb::Dict{Any,Float64}, q_ub::Dict{Any,Float64}) + + modifications = Dict{String, Any}() + modifications["per_unit"] = false + modifications["node"] = Dict{String, Any}() + modifications["pipe"] = Dict{String, Any}() + + for node in ids(wm, :node) + index = string(ref(wm, :node, node, "index")) + modifications["node"][index] = Dict("h_min" => h_lb[node], "h_max" => h_ub[node]) + end + + return modifications +end From f1750a4183634088a96d236c11b0794fa200b04d Mon Sep 17 00:00:00 2001 From: Byron Tasseff Date: Mon, 3 Aug 2020 19:55:21 -0600 Subject: [PATCH 03/24] Initial bound tightening of head. --- src/core/ref.jl | 8 +- src/util/obbt.jl | 263 ++++++++++++++++++++++++----------------------- 2 files changed, 139 insertions(+), 132 deletions(-) diff --git a/src/core/ref.jl b/src/core/ref.jl index 70cb026c..70f28ba7 100644 --- a/src/core/ref.jl +++ b/src/core/ref.jl @@ -53,13 +53,13 @@ function calc_head_bounds(wm::AbstractWaterModel, n::Int=wm.cnw) 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"]) + if haskey(node, "h_min") + head_min[i] = max(head_min[i], node["h_min"]) 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"]) + if haskey(node, "h_max") + head_max[i] = min(head_max[i], node["h_max"]) end num_junctions = length(ref(wm, n, :node_junction, i)) diff --git a/src/util/obbt.jl b/src/util/obbt.jl index 40f07d28..c0d29de1 100644 --- a/src/util/obbt.jl +++ b/src/util/obbt.jl @@ -56,6 +56,14 @@ function _get_obbt_node_var_ub(wm::AbstractWaterModel, v::Symbol) for i in ids(wm, nw, :node)) for nw in network_ids) 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 run_obbt_owf!(data::Dict{String,<:Any}, optimizer; model_type::Type = MICPRWaterModel, max_iter::Int = 100, @@ -84,11 +92,8 @@ function run_obbt_owf!(data::Dict{String,<:Any}, optimizer; Memento.error(_LOGGER, "OBBT termination criteria can only be :max or :avg.") end - # Further relax all binary variables in the model. - bin_vars = filter(v -> JuMP.is_binary(v), JuMP.all_variables(model_relaxation.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. + # Relax the model. + _relax_model!(model_relaxation) # Declare possible pass statuses. pass_statuses = [_MOI.LOCALLY_SOLVED, _MOI.OPTIMAL] @@ -119,12 +124,7 @@ function run_obbt_owf!(data::Dict{String,<:Any}, optimizer; model_bt = instantiate_model(data, model_type, WaterModels.build_mn_owf) upper_bound_constraint && _constraint_obj_bound(model_bt, upper_bound) - - # Further relax all binary variables in the bound tightening model. - bin_vars = filter(v -> JuMP.is_binary(v), JuMP.all_variables(model_bt.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. + _relax_model!(model_bt) # Initialize the statistics dictionary. stats = Dict{String,Any}() @@ -153,11 +153,14 @@ function run_obbt_owf!(data::Dict{String,<:Any}, optimizer; # Set whether or not the algorithm should immediately terminate. if termination == :avg - terminate = avg_h_reduction <= improvement_tol && avg_q_reduction <= improvement_tol + terminate = avg_h_reduction <= improvement_tol elseif termination == :max - terminate = max_h_reduction <= improvement_tol && max_q_reduction <= improvement_tol + terminate = max_h_reduction <= improvement_tol end + # Set the optimizer for model_bt. + JuMP.set_optimizer(model_bt.model, optimizer) + while !terminate # Algorithmic loop. # Set important metadata describing the current round. iter_start_time, max_h_iteration_time = time(), 0.0 @@ -165,11 +168,8 @@ function run_obbt_owf!(data::Dict{String,<:Any}, optimizer; # Loop over all subnetworks in the multinetwork. for nw in sort(collect(nw_ids(model_bt))) - # Loop over all nodes in the network. for i in ids(model_bt, nw, :node) - println(nw, " ", i) - # If the lower and upper bounds are already close, skip. if h_ub[nw][i] - h_lb[nw][i] < min_bound_width continue @@ -180,7 +180,7 @@ function run_obbt_owf!(data::Dict{String,<:Any}, optimizer; # Minimize the variable whose bounds are being tightened. JuMP.@objective(model_bt.model, _MOI.MIN_SENSE, var(model_bt, nw, :h, i)) - result_bt = optimize_model!(model_bt, optimizer=optimizer) + result_bt = optimize_model!(model_bt) # Store the new lower bound or print a warning and continue. if result_bt["termination_status"] in pass_statuses @@ -195,7 +195,7 @@ function run_obbt_owf!(data::Dict{String,<:Any}, optimizer; # Maximize the variable whose bounds are being tightened. JuMP.@objective(model_bt.model, _MOI.MAX_SENSE, var(model_bt, nw, :h, i)) - result_bt = optimize_model!(model_bt, optimizer=optimizer) + result_bt = optimize_model!(model_bt) # Store the new lower bound or print a warning and continue. if result_bt["termination_status"] in pass_statuses @@ -219,110 +219,117 @@ function run_obbt_owf!(data::Dict{String,<:Any}, optimizer; end # Revert to old bounds if new bounds are nonsensical or NaN. - (!isnan(lb) && lb > h_ub[i]) && (lb = h_lb[i]) - (!isnan(ub) && ub < h_lb[i]) && (ub = h_ub[i]) - isnan(lb) && (lb = h_lb[i]) - isnan(ub) && (ub = h_ub[i]) + (!isnan(lb) && lb > h_ub[nw][i]) && (lb = h_lb[nw][i]) + (!isnan(ub) && ub < h_lb[nw][i]) && (ub = h_ub[nw][i]) + isnan(lb) && (lb = h_lb[nw][i]) + isnan(ub) && (ub = h_ub[nw][i]) + + # Compute the h bound reduction. + h_reduction = 0.0 + + if (ub - lb >= min_bound_width) + h_reduction = (h_ub[nw][i] - h_lb[nw][i]) - (ub - lb) + h_lb[nw][i], h_ub[nw][i] = lb, ub + else + mean = 0.5 * (ub + lb) + + if mean - 0.5 * min_bound_width < h_lb[nw][i] + lb, ub = h_lb[nw][i], h_lb[nw][i] + min_bound_width + elseif mean + 0.5 * min_bound_width > h_ub[nw][i] + ub, lb = h_ub[nw][i], h_ub[nw][i] - min_bound_width + else + lb, ub = mean - 0.5 * min_bound_width, mean + 0.5 * min_bound_width + end + + h_reduction = (h_ub[nw][i] - h_lb[nw][i]) - (ub - lb) + h_lb[nw][i], h_ub[nw][i] = lb, ub + end + + total_h_reduction += h_reduction + max_h_reduction = max(h_reduction, max_h_reduction) end end + + total_count, h_range_final, avg_h_range = 0, 0.0, 0.0 + + for nw in sort(collect(nw_ids(model_bt))) + node_ids = ids(model_bt, nw, :node) + h_range_nw = sum(h_ub[nw][i] - h_lb[nw][i] for i in node_ids) + h_range_final += h_range_nw + avg_h_range += h_range_nw * inv(length(node_ids)) + total_count += length(node_ids) + end + + avg_h_reduction = total_h_reduction * inv(total_count) + parallel_time_elapsed += max_h_iteration_time + time_elapsed += time() - iter_start_time + + # Populate the modifications, update the data, and rebuild the bound tightening model. + modifications = _create_modifications(model_bt, h_lb, h_ub) + _IM.update_data!(data, modifications) + + # Further relax all binary variables in the bound tightening model. + model_bt = instantiate_model(data, model_type, WaterModels.build_mn_owf) + _relax_model!(model_bt) + + JuMP.set_optimizer(model_bt.model, optimizer) + upper_bound_constraint && _constraint_obj_bound(model_bt, upper_bound) + h = var(model_bt, :h) + + # Solve the relaxation for the updated bounds. + model_relaxation = instantiate_model(data, model_type, WaterModels.build_mn_owf) + _relax_model!(model_relaxation) + result_relaxation = optimize_model!(model_relaxation, optimizer=optimizer) + + if result_relaxation["termination_status"] in pass_statuses + current_rel_gap = (upper_bound - result_relaxation["objective"]) * inv(upper_bound) + final_relaxation_objective = result_relaxation["objective"] + else + Memento.warn(_LOGGER, "Relaxation solve failed in iteration $(current_iteration+1).") + Memento.warn(_LOGGER, "Using the previous iteration's gap to check relative gap stopping criteria.") + end + + Memento.info(_LOGGER, "Iteration $(current_iteration+1), h range: $(h_range_final), relaxation obj: $(final_relaxation_objective).") + + # Set whether or not the algorithm should terminate. + if termination == :avg + terminate = avg_h_reduction <= improvement_tol + elseif termination == :max + terminate = max_h_reduction <= improvement_tol + end + + # Iteration counter update. + current_iteration += 1 + + # Check termination criteria. + if current_iteration >= max_iter + Memento.info(_LOGGER, "Maximum iteration limit reached.") + terminate = true + elseif time_elapsed > time_limit + Memento.info(_LOGGER, "Maximum time limit reached.") + terminate = true + elseif !isinf(rel_gap_tol) && current_rel_gap < rel_gap_tol + Memento.info(_LOGGER, "Relative optimality gap < $(rel_gap_tol).") + terminate = true + end + end + + stats["final_relaxation_objective"] = final_relaxation_objective + stats["final_rel_gap_from_ub"] = isnan(upper_bound) ? Inf : current_rel_gap + stats["avg_h_range_final"], stats["h_range_final"] = 0.0, 0.0 + + for nw in sort(collect(nw_ids(model_bt))) + node_ids = ids(model_bt, nw, :node) + h_range_nw = sum(h_ub[nw][i] - h_lb[nw][i] for i in node_ids) + stats["h_range_final"] += h_range_nw + stats["avg_h_range_final"] += h_range_nw * inv(length(node_ids)) end - # # vm bound-reduction computation - # h_reduction = 0.0 - # if (ub - lb >= min_bound_width) - # h_reduction = (h_ub[node] - h_lb[node]) - (ub - lb) - # h_lb[node] = lb - # h_ub[node] = ub - # else - # mean = 0.5 * (ub + lb) - # if (mean - 0.5 * min_bound_width < h_lb[node]) - # lb = h_lb[node] - # ub = h_lb[node] + min_bound_width - # elseif (mean + 0.5 * min_bound_width > h_ub[node]) - # ub = h_ub[node] - # lb = h_ub[node] - min_bound_width - # else - # lb = mean - 0.5 * min_bound_width - # ub = mean + 0.5 * min_bound_width - # end - # h_reduction = (h_ub[node] - h_lb[node]) - (ub - lb) - # h_lb[node] = lb - # h_ub[node] = ub - # end - - # total_h_reduction += (h_reduction) - # max_h_reduction = max(h_reduction, max_h_reduction) - # end - # avg_h_reduction = total_h_reduction/length(nodes) - - # h_range_final = sum([h_ub[node] - h_lb[node] for node in nodes]) - - # q_range_final = sum([q_ub[bp] - q_lb[bp] for bp in pipes]) - - # parallel_time_elapsed += max(max_h_iteration_time, max_q_iteration_time) - - # time_elapsed += (time() - iter_start_time) - - # # populate the modifications, update the data, and rebuild the bound tightening model - # modifications = _create_modifications(model_bt, h_lb, h_ub, q_lb, q_ub) - # WaterModels.update_data!(data, modifications) - # model_bt = instantiate_model(data, model_type, WaterModels.build_mn_owf) - # (upper_bound_constraint) && (_constraint_obj_bound(model_bt, upper_bound)) - # vm = var(model_bt, :vm) - # q = var(model_bt, :q) - - # # run the qc relaxation for the updated bounds - # result_relaxation = run_owf(data, model_type::Type, optimizer) - - # if result_relaxation["termination_status"] in pass_statuses - # current_rel_gap = (upper_bound - result_relaxation["objective"])/upper_bound - # final_relaxation_objective = result_relaxation["objective"] - # else - # Memento.warn(_LOGGER, "relaxation solve failed in iteration $(current_iteration+1)") - # Memento.warn(_LOGGER, "using the previous iteration's gap to check relative gap stopping criteria") - # end - - # Memento.info(_LOGGER, "iteration $(current_iteration+1), vm range: $h_range_final, q range: $q_range_final, relaxation obj: $final_relaxation_objective") - - # # Update the termination flag. - # if termination == :avg - # terminate = avg_h_reduction <= improvement_tol && avg_q_reduction <= improvement_tol - # else - # terminate = max_h_reduction <= improvement_tol && max_q_reduction <= improvement_tol - # end - - # # iteration counter update - # current_iteration += 1 - # # check all the stopping criteria - # (current_iteration >= max_iter) && (Memento.info(_LOGGER, "maximum iteration limit reached"); break) - # (time_elapsed > time_limit) && (Memento.info(_LOGGER, "maximum time limit reached"); break) - # if (!isinf(rel_gap_tol)) && (current_rel_gap < rel_gap_tol) - # Memento.info(_LOGGER, "relative optimality gap < $rel_gap_tol") - # break - # end - - #end - - #pipes_vad_same_sign_count = 0 - #for (key, pipe) in data["pipe"] - # is_same_sign = (pipe["q_max"] >=0 && pipe["q_min"] >= 0) || (pipe["q_max"] <=0 && pipe["q_min"] <= 0) - # (is_same_sign) && (pipes_vad_same_sign_count += 1) - #end - - #stats["final_relaxation_objective"] = final_relaxation_objective - #stats["final_rel_gap_from_ub"] = isnan(upper_bound) ? Inf : current_rel_gap - #stats["h_range_final"] = h_range_final - #stats["avg_h_range_final"] = h_range_final * inv(length(nodes)) - - #stats["q_range_final"] = q_range_final - #stats["avg_q_range_final"] = q_range_final * inv(length(pipes)) - - #stats["run_time"] = time_elapsed - #stats["iteration_count"] = current_iteration - #stats["sim_parallel_run_time"] = parallel_time_elapsed - - #stats["vad_sign_determined"] = pipes_vad_same_sign_count + stats["run_time"] = time_elapsed + stats["iteration_count"] = current_iteration + stats["sim_parallel_run_time"] = parallel_time_elapsed + # Return new data dictionary and statistics. return data, stats end @@ -350,19 +357,19 @@ end function _constraint_obj_bound(wm::AbstractWaterModel, bound) end +function _create_modifications(wm::AbstractWaterModel, h_lb::Dict{Int,<:Any}, h_ub::Dict{Int,<:Any}) + modifications = Dict{String,Any}("nw"=>Dict{String,Any}(), "per_unit"=>false) -function _create_modifications(wm::AbstractWaterModel, - h_lb::Dict{Any,Float64}, h_ub::Dict{Any,Float64}, - q_lb::Dict{Any,Float64}, q_ub::Dict{Any,Float64}) + for nw in sort(collect(nw_ids(wm))) + nw_str = string(nw) + modifications["nw"][nw_str] = Dict{String,Any}() + modifications["nw"][nw_str]["node"] = Dict{String,Any}() - modifications = Dict{String, Any}() - modifications["per_unit"] = false - modifications["node"] = Dict{String, Any}() - modifications["pipe"] = Dict{String, Any}() - - for node in ids(wm, :node) - index = string(ref(wm, :node, node, "index")) - modifications["node"][index] = Dict("h_min" => h_lb[node], "h_max" => h_ub[node]) + for i in ids(wm, nw, :node) + modifications["nw"][nw_str]["node"][string(i)] = Dict{String,Any}() + modifications["nw"][nw_str]["node"][string(i)]["h_min"] = h_lb[nw][i] + modifications["nw"][nw_str]["node"][string(i)]["h_max"] = h_ub[nw][i] + end end return modifications From 394deda82fd13f32dfff4a83b15527fada7104b3 Mon Sep 17 00:00:00 2001 From: Byron Tasseff Date: Tue, 4 Aug 2020 11:48:42 -0600 Subject: [PATCH 04/24] Add node directionality (degree two) constraints. --- src/core/constraint_template.jl | 42 ++++++++++++++++++++++++++++++- src/form/directed_flow.jl | 44 ++++++++++++++++++++++++++++++++- src/form/undirected_flow.jl | 9 +++++++ src/prob/des.jl | 1 + src/prob/wf.jl | 2 ++ src/util/obbt.jl | 16 ++++++------ 6 files changed, 103 insertions(+), 11 deletions(-) diff --git a/src/core/constraint_template.jl b/src/core/constraint_template.jl index 5b0c7e95..1f97cdb4 100644 --- a/src/core/constraint_template.jl +++ b/src/core/constraint_template.jl @@ -59,6 +59,46 @@ function constraint_flow_conservation(wm::AbstractWaterModel, i::Int; nw::Int=wm 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 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) + + 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 + + ### 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`. @@ -107,7 +147,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( diff --git a/src/form/directed_flow.jl b/src/form/directed_flow.jl index 4ceceb2d..1b352157 100644 --- a/src/form/directed_flow.jl +++ b/src/form/directed_flow.jl @@ -456,6 +456,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 +528,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/undirected_flow.jl b/src/form/undirected_flow.jl index a74bfe39..689ba527 100644 --- a/src/form/undirected_flow.jl +++ b/src/form/undirected_flow.jl @@ -181,6 +181,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/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/wf.jl b/src/prob/wf.jl index af71cbc9..3e810a74 100644 --- a/src/prob/wf.jl +++ b/src/prob/wf.jl @@ -25,6 +25,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. @@ -104,6 +105,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. diff --git a/src/util/obbt.jl b/src/util/obbt.jl index c0d29de1..f2ec2697 100644 --- a/src/util/obbt.jl +++ b/src/util/obbt.jl @@ -189,7 +189,7 @@ function run_obbt_owf!(data::Dict{String,<:Any}, optimizer; lb_new = floor(10.0^precision * val) * inv(10.0^precision) lb_new > h_lb[nw][i] && (lb = lb_new) # Store the new lower bound. else - Memento.warn(_LOGGER, "BT minimization for h[$(i)] errored. Adjust tolerances.") + Memento.warn(_LOGGER, "BT minimization for $(nw)_h[$(i)] errored. Adjust tolerances.") continue end @@ -204,7 +204,7 @@ function run_obbt_owf!(data::Dict{String,<:Any}, optimizer; ub_new = ceil(10.0^precision * val) * inv(10.0^precision) ub_new < h_ub[nw][i] && (ub = ub_new) # Store the new lower bound. else - Memento.warn(_LOGGER, "BT maximization for h[$(i)] errored. Adjust tolerances.") + Memento.warn(_LOGGER, "BT maximization for $(nw)_h[$(i)] errored. Adjust tolerances.") continue end @@ -214,7 +214,7 @@ function run_obbt_owf!(data::Dict{String,<:Any}, optimizer; # Perform sanity checks on the new bounds. if lb > ub # Check if the lower bound exceeds the upper bound. - Memento.warn(_LOGGER, "BT lb > ub for h[$(i)]. Adjust tolerances.") + Memento.warn(_LOGGER, "BT lb > ub for $(nw)_h[$(i)]. Adjust tolerances.") continue end @@ -270,15 +270,14 @@ function run_obbt_owf!(data::Dict{String,<:Any}, optimizer; # Further relax all binary variables in the bound tightening model. model_bt = instantiate_model(data, model_type, WaterModels.build_mn_owf) - _relax_model!(model_bt) - + _relax_model!(model_bt) # Relax binary variables in the model. JuMP.set_optimizer(model_bt.model, optimizer) upper_bound_constraint && _constraint_obj_bound(model_bt, upper_bound) h = var(model_bt, :h) # Solve the relaxation for the updated bounds. model_relaxation = instantiate_model(data, model_type, WaterModels.build_mn_owf) - _relax_model!(model_relaxation) + _relax_model!(model_relaxation) # Relax binary variables in the model. result_relaxation = optimize_model!(model_relaxation, optimizer=optimizer) if result_relaxation["termination_status"] in pass_statuses @@ -319,10 +318,9 @@ function run_obbt_owf!(data::Dict{String,<:Any}, optimizer; stats["avg_h_range_final"], stats["h_range_final"] = 0.0, 0.0 for nw in sort(collect(nw_ids(model_bt))) - node_ids = ids(model_bt, nw, :node) - h_range_nw = sum(h_ub[nw][i] - h_lb[nw][i] for i in node_ids) + h_range_nw = sum(h_ub[nw][i] - h_lb[nw][i] for i in ids(model_bt, nw, :node)) stats["h_range_final"] += h_range_nw - stats["avg_h_range_final"] += h_range_nw * inv(length(node_ids)) + stats["avg_h_range_final"] += h_range_nw * inv(length(ids(model_bt, nw, :node))) end stats["run_time"] = time_elapsed From 5f92136ab279c1f3db3f3ad1d3a4458f6f87a465 Mon Sep 17 00:00:00 2001 From: Byron Tasseff Date: Wed, 5 Aug 2020 19:59:30 -0600 Subject: [PATCH 05/24] Fix maximum time limit termination of OBBT. --- src/util/obbt.jl | 31 +++++++++++++++++++++++-------- 1 file changed, 23 insertions(+), 8 deletions(-) diff --git a/src/util/obbt.jl b/src/util/obbt.jl index f2ec2697..005d804c 100644 --- a/src/util/obbt.jl +++ b/src/util/obbt.jl @@ -160,6 +160,7 @@ function run_obbt_owf!(data::Dict{String,<:Any}, optimizer; # Set the optimizer for model_bt. JuMP.set_optimizer(model_bt.model, optimizer) + Memento.info(_LOGGER, "Initial h range: $(stats["h_range_init"]).") while !terminate # Algorithmic loop. # Set important metadata describing the current round. @@ -170,13 +171,16 @@ function run_obbt_owf!(data::Dict{String,<:Any}, optimizer; for nw in sort(collect(nw_ids(model_bt))) # Loop over all nodes in the network. for i in ids(model_bt, nw, :node) - # If the lower and upper bounds are already close, skip. + # Capture the start time. + start_time = time() + + # If the lower and upper bounds are already sufficiently close, skip. if h_ub[nw][i] - h_lb[nw][i] < min_bound_width continue end # Store metadata for variable tightening. - lb, ub, start_time = NaN, NaN, time() + lb, ub = NaN, NaN # Minimize the variable whose bounds are being tightened. JuMP.@objective(model_bt.model, _MOI.MIN_SENSE, var(model_bt, nw, :h, i)) @@ -208,10 +212,6 @@ function run_obbt_owf!(data::Dict{String,<:Any}, optimizer; continue end - # Compute the time required for the variable's bound tightening. - end_time = time() - start_time - max_h_iteration_time = max(end_time, max_h_iteration_time) - # Perform sanity checks on the new bounds. if lb > ub # Check if the lower bound exceeds the upper bound. Memento.warn(_LOGGER, "BT lb > ub for $(nw)_h[$(i)]. Adjust tolerances.") @@ -238,7 +238,8 @@ function run_obbt_owf!(data::Dict{String,<:Any}, optimizer; elseif mean + 0.5 * min_bound_width > h_ub[nw][i] ub, lb = h_ub[nw][i], h_ub[nw][i] - min_bound_width else - lb, ub = mean - 0.5 * min_bound_width, mean + 0.5 * min_bound_width + lb = mean - 0.5 * min_bound_width + ub = mean + 0.5 * min_bound_width end h_reduction = (h_ub[nw][i] - h_lb[nw][i]) - (ub - lb) @@ -247,7 +248,22 @@ function run_obbt_owf!(data::Dict{String,<:Any}, optimizer; total_h_reduction += h_reduction max_h_reduction = max(h_reduction, max_h_reduction) + + # Compute the time required for the variable's bound tightening. + end_time = time() - start_time + max_h_iteration_time = max(end_time, max_h_iteration_time) + + # Update the time elapsed. + time_elapsed += end_time + + # Check time elapsed termination criteria. + if time_elapsed > time_limit + (terminate = true) && break + end end + + # If algorithm is supposed to terminate... + terminate && break end total_count, h_range_final, avg_h_range = 0, 0.0, 0.0 @@ -262,7 +278,6 @@ function run_obbt_owf!(data::Dict{String,<:Any}, optimizer; avg_h_reduction = total_h_reduction * inv(total_count) parallel_time_elapsed += max_h_iteration_time - time_elapsed += time() - iter_start_time # Populate the modifications, update the data, and rebuild the bound tightening model. modifications = _create_modifications(model_bt, h_lb, h_ub) From ef47db87745d8b7080740bf6b94a8c0441ee8a3d Mon Sep 17 00:00:00 2001 From: Byron Tasseff Date: Wed, 5 Aug 2020 20:37:16 -0600 Subject: [PATCH 06/24] Clean up logging of OBBT. --- src/util/obbt.jl | 77 +++++++++++++++++++++--------------------------- 1 file changed, 34 insertions(+), 43 deletions(-) diff --git a/src/util/obbt.jl b/src/util/obbt.jl index 005d804c..0f3b6261 100644 --- a/src/util/obbt.jl +++ b/src/util/obbt.jl @@ -77,8 +77,8 @@ function run_obbt_owf!(data::Dict{String,<:Any}, optimizer; termination::Symbol = :avg, kwargs...) - Memento.info(_LOGGER, "Maximum number of OBBT iterations set to default value of $(max_iter).") - Memento.info(_LOGGER, "Maximum time limit for OBBT set to default value of $(time_limit) seconds.") + Memento.info(_LOGGER, "[OBBT] Maximum number of OBBT iterations set to default value of $(max_iter).") + Memento.info(_LOGGER, "[OBBT] Maximum time limit for OBBT set to default value of $(time_limit) seconds.") model_relaxation = instantiate_model(data, model_type, WaterModels.build_mn_owf) # Check for model_type compatability with OBBT. @@ -89,7 +89,7 @@ function run_obbt_owf!(data::Dict{String,<:Any}, optimizer; # Check termination norm criteria for OBBT. if termination != :avg && termination != :max - Memento.error(_LOGGER, "OBBT termination criteria can only be :max or :avg.") + Memento.error(_LOGGER, "[OBBT] Termination criteria can only be :max or :avg.") end # Relax the model. @@ -103,15 +103,15 @@ function run_obbt_owf!(data::Dict{String,<:Any}, optimizer; current_relaxation_objective = result_relaxation["objective"] if upper_bound < current_relaxation_objective - Memento.error(_LOGGER, "The upper bound provided to OBBT is not a valid OWF upper bound.") + Memento.error(_LOGGER, "[OBBT] The provided upper bound is not valid.") end if !(result_relaxation["termination_status"] in pass_statuses) termination_status = result_relaxation["termination_status"] - Memento.warn(_LOGGER, "Initial relaxation solve status is $(termination_status).") + Memento.warn(_LOGGER, "[OBBT] Initial relaxation solve status is $(termination_status).") if termination_status == :SubOptimal - Memento.warn(_LOGGER, "Continuing with the bound tightening algorithm.") + Memento.warn(_LOGGER, "[OBBT] Continuing with the bound tightening algorithm.") end end @@ -119,7 +119,7 @@ function run_obbt_owf!(data::Dict{String,<:Any}, optimizer; if !isinf(upper_bound) current_rel_gap = (upper_bound - current_relaxation_objective) * inv(upper_bound) - Memento.info(_LOGGER, "Initial relaxation gap is $(current_rel_gap).") + Memento.info(_LOGGER, "[OBBT] Initial relaxation gap is $(current_rel_gap).") end model_bt = instantiate_model(data, model_type, WaterModels.build_mn_owf) @@ -135,41 +135,37 @@ function run_obbt_owf!(data::Dict{String,<:Any}, optimizer; h_lb = _get_obbt_node_var_lb(model_bt, :h) h_ub = _get_obbt_node_var_ub(model_bt, :h) - stats["h_range_init"], stats["avg_h_range_init"] = 0.0, 0.0 + total_count, h_range, avg_h_range = 0, 0.0, 0.0 for nw in sort(collect(nw_ids(model_bt))) node_ids = ids(model_bt, nw, :node) h_range_nw = sum(h_ub[nw][i] - h_lb[nw][i] for i in node_ids) - stats["h_range_init"] += h_range_nw - stats["avg_h_range_init"] += h_range_nw * inv(length(node_ids)) + h_range += h_range_nw + total_count += length(node_ids) end + avg_h_range = h_range * inv(total_count) + stats["h_range_init"] = h_range + stats["avg_h_range_init"] = avg_h_range + # Initialize algorithm termination metadata. - h_range_final, total_h_reduction = 0.0, Inf - max_h_reduction, avg_h_reduction = Inf, Inf - final_relaxation_objective = NaN current_iteration, time_elapsed = 0, 0.0 parallel_time_elapsed, terminate = 0.0, false - - # Set whether or not the algorithm should immediately terminate. - if termination == :avg - terminate = avg_h_reduction <= improvement_tol - elseif termination == :max - terminate = max_h_reduction <= improvement_tol - end + final_relaxation_objective = NaN # Set the optimizer for model_bt. JuMP.set_optimizer(model_bt.model, optimizer) - Memento.info(_LOGGER, "Initial h range: $(stats["h_range_init"]).") + Memento.info(_LOGGER, "[OBBT] Initial average h range: $(stats["avg_h_range_init"]).") while !terminate # Algorithmic loop. # Set important metadata describing the current round. iter_start_time, max_h_iteration_time = time(), 0.0 - total_h_reduction, avg_h_reduction, max_h_reduction = 0.0, 0.0, 0.0 + total_h_reduction, max_h_reduction = 0.0, 0.0 # Loop over all subnetworks in the multinetwork. for nw in sort(collect(nw_ids(model_bt))) # Loop over all nodes in the network. + for i in ids(model_bt, nw, :node) # Capture the start time. start_time = time() @@ -193,7 +189,7 @@ function run_obbt_owf!(data::Dict{String,<:Any}, optimizer; lb_new = floor(10.0^precision * val) * inv(10.0^precision) lb_new > h_lb[nw][i] && (lb = lb_new) # Store the new lower bound. else - Memento.warn(_LOGGER, "BT minimization for $(nw)_h[$(i)] errored. Adjust tolerances.") + Memento.warn(_LOGGER, "[OBBT] Minimization for $(nw)_h[$(i)] errored. Adjust tolerances.") continue end @@ -208,13 +204,13 @@ function run_obbt_owf!(data::Dict{String,<:Any}, optimizer; ub_new = ceil(10.0^precision * val) * inv(10.0^precision) ub_new < h_ub[nw][i] && (ub = ub_new) # Store the new lower bound. else - Memento.warn(_LOGGER, "BT maximization for $(nw)_h[$(i)] errored. Adjust tolerances.") + Memento.warn(_LOGGER, "[OBBT] Maximization for $(nw)_h[$(i)] errored. Adjust tolerances.") continue end # Perform sanity checks on the new bounds. if lb > ub # Check if the lower bound exceeds the upper bound. - Memento.warn(_LOGGER, "BT lb > ub for $(nw)_h[$(i)]. Adjust tolerances.") + Memento.warn(_LOGGER, "[OBBT] lb > ub for $(nw)_h[$(i)]. Adjust tolerances.") continue end @@ -266,16 +262,16 @@ function run_obbt_owf!(data::Dict{String,<:Any}, optimizer; terminate && break end - total_count, h_range_final, avg_h_range = 0, 0.0, 0.0 + total_count, h_range = 0, 0.0 for nw in sort(collect(nw_ids(model_bt))) node_ids = ids(model_bt, nw, :node) h_range_nw = sum(h_ub[nw][i] - h_lb[nw][i] for i in node_ids) - h_range_final += h_range_nw - avg_h_range += h_range_nw * inv(length(node_ids)) + h_range += h_range_nw total_count += length(node_ids) end + avg_h_range = h_range * inv(total_count) avg_h_reduction = total_h_reduction * inv(total_count) parallel_time_elapsed += max_h_iteration_time @@ -299,11 +295,11 @@ function run_obbt_owf!(data::Dict{String,<:Any}, optimizer; current_rel_gap = (upper_bound - result_relaxation["objective"]) * inv(upper_bound) final_relaxation_objective = result_relaxation["objective"] else - Memento.warn(_LOGGER, "Relaxation solve failed in iteration $(current_iteration+1).") - Memento.warn(_LOGGER, "Using the previous iteration's gap to check relative gap stopping criteria.") + Memento.warn(_LOGGER, "[OBBT] Relaxation solve failed in iteration $(current_iteration+1).") + Memento.warn(_LOGGER, "[OBBT] Using the previous iteration's gap to check relative gap stopping criteria.") end - Memento.info(_LOGGER, "Iteration $(current_iteration+1), h range: $(h_range_final), relaxation obj: $(final_relaxation_objective).") + Memento.info(_LOGGER, "[OBBT] Iteration $(current_iteration+1): average h range: $(avg_h_range).") # Set whether or not the algorithm should terminate. if termination == :avg @@ -317,26 +313,21 @@ function run_obbt_owf!(data::Dict{String,<:Any}, optimizer; # Check termination criteria. if current_iteration >= max_iter - Memento.info(_LOGGER, "Maximum iteration limit reached.") + Memento.info(_LOGGER, "[OBBT] Maximum iteration limit reached.") terminate = true elseif time_elapsed > time_limit - Memento.info(_LOGGER, "Maximum time limit reached.") + Memento.info(_LOGGER, "[OBBT] Maximum time limit reached.") terminate = true elseif !isinf(rel_gap_tol) && current_rel_gap < rel_gap_tol - Memento.info(_LOGGER, "Relative optimality gap < $(rel_gap_tol).") + Memento.info(_LOGGER, "[OBBT] Relative optimality gap < $(rel_gap_tol).") terminate = true end end + stats["h_range_final"] = h_range + stats["avg_h_range_final"] = avg_h_range stats["final_relaxation_objective"] = final_relaxation_objective stats["final_rel_gap_from_ub"] = isnan(upper_bound) ? Inf : current_rel_gap - stats["avg_h_range_final"], stats["h_range_final"] = 0.0, 0.0 - - for nw in sort(collect(nw_ids(model_bt))) - h_range_nw = sum(h_ub[nw][i] - h_lb[nw][i] for i in ids(model_bt, nw, :node)) - stats["h_range_final"] += h_range_nw - stats["avg_h_range_final"] += h_range_nw * inv(length(ids(model_bt, nw, :node))) - end stats["run_time"] = time_elapsed stats["iteration_count"] = current_iteration @@ -358,11 +349,11 @@ end function _check_obbt_options(ub::Float64, rel_gap::Float64, ub_constraint::Bool) if ub_constraint && isinf(ub) - Memento.error(_LOGGER, "The option \"upper_bound_constraint\" cannot be set to true without specifying an upper bound.") + Memento.error(_LOGGER, "[OBBT] The option \"upper_bound_constraint\" cannot be set to true without specifying an upper bound.") end if !isinf(rel_gap) && isinf(ub) - Memento.error(_LOGGER, "The option \"rel_gap_tol\" is specified without providing an upper bound.") + Memento.error(_LOGGER, "[OBBT] The option \"rel_gap_tol\" is specified without providing an upper bound.") end end From 64eb092e8292e66ee3c12186e95e09ff65367312 Mon Sep 17 00:00:00 2001 From: Byron Tasseff Date: Wed, 5 Aug 2020 22:04:48 -0600 Subject: [PATCH 07/24] Add more OBBT options. --- src/util/obbt.jl | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/src/util/obbt.jl b/src/util/obbt.jl index 0f3b6261..da2aaf02 100644 --- a/src/util/obbt.jl +++ b/src/util/obbt.jl @@ -75,11 +75,13 @@ function run_obbt_owf!(data::Dict{String,<:Any}, optimizer; improvement_tol::Float64 = 1.0e-3, precision::Int = 4, termination::Symbol = :avg, + relaxed::Bool = true, + ext::Dict{Symbol,<:Any} = Dict{Symbol,Any}(:pump_breakpoints=>5), kwargs...) Memento.info(_LOGGER, "[OBBT] Maximum number of OBBT iterations set to default value of $(max_iter).") Memento.info(_LOGGER, "[OBBT] Maximum time limit for OBBT set to default value of $(time_limit) seconds.") - model_relaxation = instantiate_model(data, model_type, WaterModels.build_mn_owf) + model_relaxation = instantiate_model(data, model_type, WaterModels.build_mn_owf; ext=ext) # Check for model_type compatability with OBBT. _check_variables(model_relaxation) @@ -122,9 +124,9 @@ function run_obbt_owf!(data::Dict{String,<:Any}, optimizer; Memento.info(_LOGGER, "[OBBT] Initial relaxation gap is $(current_rel_gap).") end - model_bt = instantiate_model(data, model_type, WaterModels.build_mn_owf) + model_bt = instantiate_model(data, model_type, WaterModels.build_mn_owf; ext=ext) upper_bound_constraint && _constraint_obj_bound(model_bt, upper_bound) - _relax_model!(model_bt) + relaxed && _relax_model!(model_bt) # Initialize the statistics dictionary. stats = Dict{String,Any}() @@ -280,15 +282,15 @@ function run_obbt_owf!(data::Dict{String,<:Any}, optimizer; _IM.update_data!(data, modifications) # Further relax all binary variables in the bound tightening model. - model_bt = instantiate_model(data, model_type, WaterModels.build_mn_owf) - _relax_model!(model_bt) # Relax binary variables in the model. + model_bt = instantiate_model(data, model_type, WaterModels.build_mn_owf; ext=ext) + relaxed && _relax_model!(model_bt) # Relax binary variables in the model. JuMP.set_optimizer(model_bt.model, optimizer) upper_bound_constraint && _constraint_obj_bound(model_bt, upper_bound) h = var(model_bt, :h) # Solve the relaxation for the updated bounds. - model_relaxation = instantiate_model(data, model_type, WaterModels.build_mn_owf) - _relax_model!(model_relaxation) # Relax binary variables in the model. + model_relaxation = instantiate_model(data, model_type, WaterModels.build_mn_owf; ext=ext) + relaxed && _relax_model!(model_relaxation) # Relax binary variables in the model. result_relaxation = optimize_model!(model_relaxation, optimizer=optimizer) if result_relaxation["termination_status"] in pass_statuses From 30b0ba8bd8af365ea25c159991616b54a78282b6 Mon Sep 17 00:00:00 2001 From: Byron Tasseff Date: Mon, 10 Aug 2020 14:51:29 -0600 Subject: [PATCH 08/24] Simplifying OBBT utility. --- src/core/ref.jl | 28 +-- src/form/directed_flow.jl | 9 - src/util/obbt.jl | 433 ++++++++++++++------------------------ 3 files changed, 173 insertions(+), 297 deletions(-) diff --git a/src/core/ref.jl b/src/core/ref.jl index 70f28ba7..efb9e22d 100644 --- a/src/core/ref.jl +++ b/src/core/ref.jl @@ -52,16 +52,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, "h_min") - head_min[i] = max(head_min[i], node["h_min"]) - end - - # Override the maximum head value at node i, if additional data is present. - if haskey(node, "h_max") - head_max[i] = min(head_max[i], node["h_max"]) - 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)) @@ -93,6 +83,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 @@ -144,6 +139,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 +149,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 +162,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/form/directed_flow.jl b/src/form/directed_flow.jl index 1b352157..51977ed5 100644 --- a/src/form/directed_flow.jl +++ b/src/form/directed_flow.jl @@ -190,9 +190,6 @@ 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. c_1 = JuMP.@constraint(wm.model, qp <= JuMP.upper_bound(qp) * z) c_2 = JuMP.@constraint(wm.model, qp >= _q_eps * z) @@ -286,9 +283,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) @@ -340,9 +334,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) diff --git a/src/util/obbt.jl b/src/util/obbt.jl index da2aaf02..a5e10e21 100644 --- a/src/util/obbt.jl +++ b/src/util/obbt.jl @@ -32,7 +32,7 @@ data, stats = run_obbt_owf!("examples/data/epanet/van_zyl.inp", ipopt) * `rel_gap_tol`: tolerance used to terminate the algorithm when the objective value of the relaxation is close to the upper bound specified using the `upper_bound` keyword. -* `min_bound_width`: domain beyond which bound tightening is not performed. +* `min_width`: domain beyond which bound tightening is not performed. * `termination`: Bound-tightening algorithm terminates if the improvement in the average or maximum bound improvement, specified using either the `termination = :avg` or the `termination = :max` option, is less than @@ -44,17 +44,6 @@ function run_obbt_owf!(file::String, optimizer; kwargs...) return run_obbt_owf!(data, optimizer; kwargs...) end -function _get_obbt_node_var_lb(wm::AbstractWaterModel, v::Symbol) - network_ids = sort(collect(nw_ids(wm))) - return Dict(nw => Dict(i => JuMP.lower_bound(var(wm, nw, v, i)) - for i in ids(wm, nw, :node)) for nw in network_ids) -end - -function _get_obbt_node_var_ub(wm::AbstractWaterModel, v::Symbol) - network_ids = sort(collect(nw_ids(wm))) - return Dict(nw => Dict(i => JuMP.upper_bound(var(wm, nw, v, i)) - for i in ids(wm, nw, :node)) for nw in network_ids) -end function _relax_model!(wm::AbstractWaterModel) # Further relax all binary variables in the model. @@ -64,287 +53,160 @@ function _relax_model!(wm::AbstractWaterModel) JuMP.set_upper_bound.(bin_vars, 1.0) # Upper-bound the relaxed binary variables. end -function run_obbt_owf!(data::Dict{String,<:Any}, optimizer; - model_type::Type = MICPRWaterModel, - max_iter::Int = 100, - time_limit::Float64 = 3600.0, - upper_bound::Float64 = Inf, - upper_bound_constraint::Bool = false, - rel_gap_tol::Float64 = Inf, - min_bound_width::Float64 = 1.0e-2, - improvement_tol::Float64 = 1.0e-3, - precision::Int = 4, - termination::Symbol = :avg, - relaxed::Bool = true, - ext::Dict{Symbol,<:Any} = Dict{Symbol,Any}(:pump_breakpoints=>5), - kwargs...) - - Memento.info(_LOGGER, "[OBBT] Maximum number of OBBT iterations set to default value of $(max_iter).") - Memento.info(_LOGGER, "[OBBT] Maximum time limit for OBBT set to default value of $(time_limit) seconds.") - model_relaxation = instantiate_model(data, model_type, WaterModels.build_mn_owf; ext=ext) - # Check for model_type compatability with OBBT. - _check_variables(model_relaxation) +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 - # Check for other keyword argument consistencies. - _check_obbt_options(upper_bound, rel_gap_tol, upper_bound_constraint) - # Check termination norm criteria for OBBT. - if termination != :avg && termination != :max - Memento.error(_LOGGER, "[OBBT] Termination criteria can only be :max or :avg.") - 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 - # Relax the model. - _relax_model!(model_relaxation) - # Declare possible pass statuses. - pass_statuses = [_MOI.LOCALLY_SOLVED, _MOI.OPTIMAL] +function _create_modifications(wm::AbstractWaterModel) + data = Dict{String,Any}("nw"=>Dict{String,Any}(), "per_unit"=>false) - # Compute initial relative gap between relaxation objective and upper_bound. - result_relaxation = optimize_model!(model_relaxation, optimizer=optimizer) - current_relaxation_objective = result_relaxation["objective"] + 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 - if upper_bound < current_relaxation_objective - Memento.error(_LOGGER, "[OBBT] The provided upper bound is not valid.") + 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 - if !(result_relaxation["termination_status"] in pass_statuses) - termination_status = result_relaxation["termination_status"] - Memento.warn(_LOGGER, "[OBBT] Initial relaxation solve status is $(termination_status).") + return data +end - if termination_status == :SubOptimal - Memento.warn(_LOGGER, "[OBBT] Continuing with the bound tightening algorithm.") - end - end - current_rel_gap = Inf +function _get_head_index_set(wm::AbstractWaterModel; width::Float64=0.1, prec::Int=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-5, prec::Int=7) + 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]) + data_index = data["nw"][nw_str][comp_str][index_str] + return data_index[var_str * "_min"], data_index[var_str * "_max"] +end - if !isinf(upper_bound) - current_rel_gap = (upper_bound - current_relaxation_objective) * inv(upper_bound) - Memento.info(_LOGGER, "[OBBT] Initial relaxation gap is $(current_rel_gap).") + +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 + 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] end +end + + +function _print_average_widths(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)), " - model_bt = instantiate_model(data, model_type, WaterModels.build_mn_owf; ext=ext) - upper_bound_constraint && _constraint_obj_bound(model_bt, upper_bound) - relaxed && _relax_model!(model_bt) - - # Initialize the statistics dictionary. - stats = Dict{String,Any}() - stats["model_type"] = model_type - stats["initial_relaxation_objective"] = current_relaxation_objective - stats["initial_rel_gap_from_ub"] = current_rel_gap - stats["upper_bound"] = upper_bound - - h_lb = _get_obbt_node_var_lb(model_bt, :h) - h_ub = _get_obbt_node_var_ub(model_bt, :h) - total_count, h_range, avg_h_range = 0, 0.0, 0.0 - - for nw in sort(collect(nw_ids(model_bt))) - node_ids = ids(model_bt, nw, :node) - h_range_nw = sum(h_ub[nw][i] - h_lb[nw][i] for i in node_ids) - h_range += h_range_nw - total_count += length(node_ids) + 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)), " + end end - avg_h_range = h_range * inv(total_count) - stats["h_range_init"] = h_range - stats["avg_h_range_init"] = avg_h_range + Memento.info(_LOGGER, message[1:end-2] * ".") +end - # Initialize algorithm termination metadata. - current_iteration, time_elapsed = 0, 0.0 - parallel_time_elapsed, terminate = 0.0, false - final_relaxation_objective = NaN - # Set the optimizer for model_bt. - JuMP.set_optimizer(model_bt.model, optimizer) - Memento.info(_LOGGER, "[OBBT] Initial average h range: $(stats["avg_h_range_init"]).") +function run_obbt_owf!(data::Dict{String,<:Any}, optimizer; + model_type::Type = MICPRWaterModel, time_limit::Float64 = 3600.0, upper_bound::Float64 = + Inf, upper_bound_constraint::Bool = false, rel_gap_tol::Float64 = Inf, + improvement_tol::Float64 = 1.0e-3, termination::Symbol = :avg, relaxed::Bool = true, + 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.") - while !terminate # Algorithmic loop. - # Set important metadata describing the current round. - iter_start_time, max_h_iteration_time = time(), 0.0 - total_h_reduction, max_h_reduction = 0.0, 0.0 - - # Loop over all subnetworks in the multinetwork. - for nw in sort(collect(nw_ids(model_bt))) - # Loop over all nodes in the network. - - for i in ids(model_bt, nw, :node) - # Capture the start time. - start_time = time() - - # If the lower and upper bounds are already sufficiently close, skip. - if h_ub[nw][i] - h_lb[nw][i] < min_bound_width - continue - end - - # Store metadata for variable tightening. - lb, ub = NaN, NaN - - # Minimize the variable whose bounds are being tightened. - JuMP.@objective(model_bt.model, _MOI.MIN_SENSE, var(model_bt, nw, :h, i)) - result_bt = optimize_model!(model_bt) - - # Store the new lower bound or print a warning and continue. - if result_bt["termination_status"] in pass_statuses - # Compute the candidate lower bound at the specified precision. - val = JuMP.objective_value(model_bt.model) - lb_new = floor(10.0^precision * val) * inv(10.0^precision) - lb_new > h_lb[nw][i] && (lb = lb_new) # Store the new lower bound. - else - Memento.warn(_LOGGER, "[OBBT] Minimization for $(nw)_h[$(i)] errored. Adjust tolerances.") - continue - end - - # Maximize the variable whose bounds are being tightened. - JuMP.@objective(model_bt.model, _MOI.MAX_SENSE, var(model_bt, nw, :h, i)) - result_bt = optimize_model!(model_bt) - - # Store the new lower bound or print a warning and continue. - if result_bt["termination_status"] in pass_statuses - # Compute the candidate lower bound at the specified precision. - val = JuMP.objective_value(model_bt.model) - ub_new = ceil(10.0^precision * val) * inv(10.0^precision) - ub_new < h_ub[nw][i] && (ub = ub_new) # Store the new lower bound. - else - Memento.warn(_LOGGER, "[OBBT] Maximization for $(nw)_h[$(i)] errored. Adjust tolerances.") - continue - end - - # Perform sanity checks on the new bounds. - if lb > ub # Check if the lower bound exceeds the upper bound. - Memento.warn(_LOGGER, "[OBBT] lb > ub for $(nw)_h[$(i)]. Adjust tolerances.") - continue - end - - # Revert to old bounds if new bounds are nonsensical or NaN. - (!isnan(lb) && lb > h_ub[nw][i]) && (lb = h_lb[nw][i]) - (!isnan(ub) && ub < h_lb[nw][i]) && (ub = h_ub[nw][i]) - isnan(lb) && (lb = h_lb[nw][i]) - isnan(ub) && (ub = h_ub[nw][i]) - - # Compute the h bound reduction. - h_reduction = 0.0 - - if (ub - lb >= min_bound_width) - h_reduction = (h_ub[nw][i] - h_lb[nw][i]) - (ub - lb) - h_lb[nw][i], h_ub[nw][i] = lb, ub - else - mean = 0.5 * (ub + lb) - - if mean - 0.5 * min_bound_width < h_lb[nw][i] - lb, ub = h_lb[nw][i], h_lb[nw][i] + min_bound_width - elseif mean + 0.5 * min_bound_width > h_ub[nw][i] - ub, lb = h_ub[nw][i], h_ub[nw][i] - min_bound_width - else - lb = mean - 0.5 * min_bound_width - ub = mean + 0.5 * min_bound_width - end - - h_reduction = (h_ub[nw][i] - h_lb[nw][i]) - (ub - lb) - h_lb[nw][i], h_ub[nw][i] = lb, ub - end - - total_h_reduction += h_reduction - max_h_reduction = max(h_reduction, max_h_reduction) - - # Compute the time required for the variable's bound tightening. - end_time = time() - start_time - max_h_iteration_time = max(end_time, max_h_iteration_time) - - # Update the time elapsed. - time_elapsed += end_time - - # Check time elapsed termination criteria. - if time_elapsed > time_limit - (terminate = true) && break - end - end - - # If algorithm is supposed to terminate... - terminate && break - end + # Check for keyword argument inconsistencies. + _check_obbt_options(upper_bound, rel_gap_tol, upper_bound_constraint) - total_count, h_range = 0, 0.0 + # Instantiate the bound tightening model and relax integrality, if specified. + bt = instantiate_model(data, model_type, WaterModels.build_mn_owf; ext=ext) + upper_bound_constraint && _constraint_obj_bound(bt, upper_bound) + relaxed && _relax_model!(bt) - for nw in sort(collect(nw_ids(model_bt))) - node_ids = ids(model_bt, nw, :node) - h_range_nw = sum(h_ub[nw][i] - h_lb[nw][i] for i in node_ids) - h_range += h_range_nw - total_count += length(node_ids) - end + # Build the dictionary and sets that will store to the network. + modifications = _create_modifications(bt) + var_index_set = vcat(_get_head_index_set(bt), _get_flow_index_set(bt)) + bounds = [_get_existing_bounds(modifications, vid) for vid in var_index_set] - avg_h_range = h_range * inv(total_count) - avg_h_reduction = total_h_reduction * inv(total_count) - parallel_time_elapsed += max_h_iteration_time + # Print the initial average bound widths. + _print_average_widths(modifications) - # Populate the modifications, update the data, and rebuild the bound tightening model. - modifications = _create_modifications(model_bt, h_lb, h_ub) - _IM.update_data!(data, modifications) + # Initialize algorithm termination metadata. + current_iteration, time_elapsed = 1, 0.0 + parallel_time_elapsed, terminate = 0.0, false - # Further relax all binary variables in the bound tightening model. - model_bt = instantiate_model(data, model_type, WaterModels.build_mn_owf; ext=ext) - relaxed && _relax_model!(model_bt) # Relax binary variables in the model. - JuMP.set_optimizer(model_bt.model, optimizer) - upper_bound_constraint && _constraint_obj_bound(model_bt, upper_bound) - h = var(model_bt, :h) - - # Solve the relaxation for the updated bounds. - model_relaxation = instantiate_model(data, model_type, WaterModels.build_mn_owf; ext=ext) - relaxed && _relax_model!(model_relaxation) # Relax binary variables in the model. - result_relaxation = optimize_model!(model_relaxation, optimizer=optimizer) - - if result_relaxation["termination_status"] in pass_statuses - current_rel_gap = (upper_bound - result_relaxation["objective"]) * inv(upper_bound) - final_relaxation_objective = result_relaxation["objective"] - else - Memento.warn(_LOGGER, "[OBBT] Relaxation solve failed in iteration $(current_iteration+1).") - Memento.warn(_LOGGER, "[OBBT] Using the previous iteration's gap to check relative gap stopping criteria.") - end + # Set the optimizer for the bound tightening model. + JuMP.set_optimizer(bt.model, optimizer) - Memento.info(_LOGGER, "[OBBT] Iteration $(current_iteration+1): average h range: $(avg_h_range).") + while !terminate # Algorithmic loop. + Memento.info(_LOGGER, "[OBBT] Starting iteration $(current_iteration).") - # Set whether or not the algorithm should terminate. - if termination == :avg - terminate = avg_h_reduction <= improvement_tol - elseif termination == :max - terminate = max_h_reduction <= improvement_tol - end + # 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 - # Iteration counter update. - current_iteration += 1 + # 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) - # Check termination criteria. - if current_iteration >= max_iter - Memento.info(_LOGGER, "[OBBT] Maximum iteration limit reached.") - terminate = true - elseif time_elapsed > time_limit - Memento.info(_LOGGER, "[OBBT] Maximum time limit reached.") - terminate = true - elseif !isinf(rel_gap_tol) && current_rel_gap < rel_gap_tol - Memento.info(_LOGGER, "[OBBT] Relative optimality gap < $(rel_gap_tol).") - terminate = true + # 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 - end - - stats["h_range_final"] = h_range - stats["avg_h_range_final"] = avg_h_range - stats["final_relaxation_objective"] = final_relaxation_objective - stats["final_rel_gap_from_ub"] = isnan(upper_bound) ? Inf : current_rel_gap - stats["run_time"] = time_elapsed - stats["iteration_count"] = current_iteration - stats["sim_parallel_run_time"] = parallel_time_elapsed + # Update the iteration counter. + current_iteration += 1 - # Return new data dictionary and statistics. - return data, stats -end + # 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) + # Print the current average bound widths. + _print_average_widths(data) -function _check_variables(wm::AbstractWaterModel) - try - h = var(wm, :h) - catch err - (isa(error, KeyError)) && (Memento.error(_LOGGER, "OBBT is not supported for models without head variables.")) + # Instantiate the new model and set the optimizer. + bt = instantiate_model(data, model_type, build_mn_owf; ext=ext) + upper_bound_constraint && _constraint_obj_bound(bt, upper_bound) + relaxed && _relax_model!(bt) + JuMP.set_optimizer(bt.model, optimizer) end end @@ -360,23 +222,44 @@ function _check_obbt_options(ub::Float64, rel_gap::Float64, ub_constraint::Bool) end -function _constraint_obj_bound(wm::AbstractWaterModel, bound) -end +function _compute_corrected_bounds(lb::Float64, ub::Float64, lb_old::Float64, ub_old::Float64, width::Float64, prec::Int) + # Convert new bounds to the desired decimal precision. + lb = max(floor(10.0^prec * lb) * inv(10.0^prec), lb_old) + ub = min(ceil(10.0^prec * ub) * inv(10.0^prec), ub_old) -function _create_modifications(wm::AbstractWaterModel, h_lb::Dict{Int,<:Any}, h_ub::Dict{Int,<:Any}) - modifications = Dict{String,Any}("nw"=>Dict{String,Any}(), "per_unit"=>false) + # If the bounds satisfy the minimum width, return them. + ub - lb >= width && return lb, ub - for nw in sort(collect(nw_ids(wm))) - nw_str = string(nw) - modifications["nw"][nw_str] = Dict{String,Any}() - modifications["nw"][nw_str]["node"] = Dict{String,Any}() - - for i in ids(wm, nw, :node) - modifications["nw"][nw_str]["node"][string(i)] = Dict{String,Any}() - modifications["nw"][nw_str]["node"][string(i)]["h_min"] = h_lb[nw][i] - modifications["nw"][nw_str]["node"][string(i)]["h_max"] = h_ub[nw][i] - end + # 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 + - return modifications +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 $(JuMP.name(variable)) 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 From a3a2999d58fa492349a78abbca2ff8fb4304c76f Mon Sep 17 00:00:00 2001 From: Byron Tasseff Date: Mon, 10 Aug 2020 15:50:22 -0600 Subject: [PATCH 09/24] Fix logging error. --- src/util/obbt.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/util/obbt.jl b/src/util/obbt.jl index a5e10e21..ec7beda2 100644 --- a/src/util/obbt.jl +++ b/src/util/obbt.jl @@ -258,7 +258,7 @@ function _solve_bound(wm::AbstractWaterModel, index::Tuple, sense::_MOI.Optimiza candidate = JuMP.objective_value(wm.model) return sense === _MOI.MIN_SENSE ? max(candidate, start) : min(candidate, start) else - message = "[OBBT] Optimization of $(JuMP.name(variable)) errored. Adjust tolerances." + 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 From be2d5e7e9cf9f89658e1082e7901a64a2cbb32f3 Mon Sep 17 00:00:00 2001 From: Byron Tasseff Date: Mon, 10 Aug 2020 15:56:50 -0600 Subject: [PATCH 10/24] Fix logging error. --- src/util/obbt.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/util/obbt.jl b/src/util/obbt.jl index ec7beda2..379da036 100644 --- a/src/util/obbt.jl +++ b/src/util/obbt.jl @@ -258,7 +258,7 @@ function _solve_bound(wm::AbstractWaterModel, index::Tuple, sense::_MOI.Optimiza 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." + 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 From 80f05a124d01874658c96e0942bd90b69798af55 Mon Sep 17 00:00:00 2001 From: Byron Tasseff Date: Tue, 11 Aug 2020 12:32:57 -0600 Subject: [PATCH 11/24] Remove epsilon requirement from check valve formulation. --- src/form/directed_flow.jl | 17 ++++++++--------- src/form/undirected_flow.jl | 9 ++++----- test/wf.jl | 1 - 3 files changed, 12 insertions(+), 15 deletions(-) diff --git a/src/form/directed_flow.jl b/src/form/directed_flow.jl index 51977ed5..d34d7f01 100644 --- a/src/form/directed_flow.jl +++ b/src/form/directed_flow.jl @@ -190,9 +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) - # 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) @@ -200,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 diff --git a/src/form/undirected_flow.jl b/src/form/undirected_flow.jl index 689ba527..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) * min(0.0, 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 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 From 1723a6889223dfec668dc7a410aaf7f9942f4ec5 Mon Sep 17 00:00:00 2001 From: Byron Tasseff Date: Tue, 11 Aug 2020 22:09:31 -0600 Subject: [PATCH 12/24] Fix bug in reservoir head time series. --- src/core/constraint_template.jl | 3 ++- src/core/objective.jl | 9 --------- src/core/ref.jl | 2 +- src/form/milpr.jl | 25 +++++++++++++++---------- src/io/epanet.jl | 2 +- 5 files changed, 19 insertions(+), 22 deletions(-) diff --git a/src/core/constraint_template.jl b/src/core/constraint_template.jl index 1f97cdb4..065f96ae 100644 --- a/src/core/constraint_template.jl +++ b/src/core/constraint_template.jl @@ -78,10 +78,11 @@ function constraint_node_directionality(wm::AbstractWaterModel, i::Int; nw::Int= reservoirs = ref(wm, nw, :node_reservoir) num_components = length(junctions) + length(tanks) + length(reservoirs) - # Get the degree of node `i`. + # 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) 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 efb9e22d..cbdababa 100644 --- a/src/core/ref.jl +++ b/src/core/ref.jl @@ -98,7 +98,7 @@ function calc_flow_bounds(wm::AbstractWaterModel, n::Int=wm.cnw) alpha = ref(wm, n, :alpha) junctions = values(ref(wm, n, :junction)) - sum_demand = sum(abs(junction["demand"]) for junction in junctions) + sum_demand = sum(junction["demand"] for junction in junctions) if :time_step in keys(ref(wm, n)) time_step = ref(wm, n, :time_step) diff --git a/src/form/milpr.jl b/src/form/milpr.jl index 963ff302..395edbd7 100644 --- a/src/form/milpr.jl +++ b/src/form/milpr.jl @@ -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_oa_z(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 @@ -158,7 +163,7 @@ function constraint_shutoff_valve_head_loss(wm::AbstractMILPRModel, n::Int, a::I # 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)) + 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) diff --git a/src/io/epanet.jl b/src/io/epanet.jl index 5da33a6a..a8293010 100644 --- a/src/io/epanet.jl +++ b/src/io/epanet.jl @@ -282,7 +282,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. From 4968123bd37807fc2d741780de97267bf350fb88 Mon Sep 17 00:00:00 2001 From: Byron Tasseff Date: Tue, 18 Aug 2020 11:26:56 -0600 Subject: [PATCH 13/24] Beginning work on #100, #101. --- src/core/constraint_template.jl | 38 ++++++++++++++++++++------------- src/core/ref.jl | 10 +++++++-- src/io/epanet.jl | 10 ++++++++- src/prob/wf.jl | 8 +++---- src/util/obbt.jl | 5 ++++- 5 files changed, 48 insertions(+), 23 deletions(-) diff --git a/src/core/constraint_template.jl b/src/core/constraint_template.jl index 065f96ae..90e0acb6 100644 --- a/src/core/constraint_template.jl +++ b/src/core/constraint_template.jl @@ -127,10 +127,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 @@ -164,13 +167,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 @@ -179,10 +184,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_1, :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/ref.jl b/src/core/ref.jl index cbdababa..0d7893ae 100644 --- a/src/core/ref.jl +++ b/src/core/ref.jl @@ -65,8 +65,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) diff --git a/src/io/epanet.jl b/src/io/epanet.jl index a8293010..89b6ec31 100644 --- a/src/io/epanet.jl +++ b/src/io/epanet.jl @@ -668,6 +668,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" @@ -991,6 +992,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" @@ -1046,6 +1048,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" @@ -1123,11 +1126,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 @@ -1145,6 +1150,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/wf.jl b/src/prob/wf.jl index 3e810a74..5b9f6f47 100644 --- a/src/prob/wf.jl +++ b/src/prob/wf.jl @@ -61,9 +61,9 @@ 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"] < 0.0 constraint_source_directionality(wm, junction["node"]) end end @@ -141,9 +141,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 index 379da036..90319866 100644 --- a/src/util/obbt.jl +++ b/src/util/obbt.jl @@ -145,7 +145,7 @@ end function run_obbt_owf!(data::Dict{String,<:Any}, optimizer; model_type::Type = MICPRWaterModel, time_limit::Float64 = 3600.0, upper_bound::Float64 = - Inf, upper_bound_constraint::Bool = false, rel_gap_tol::Float64 = Inf, + Inf, upper_bound_constraint::Bool = false, rel_gap_tol::Float64 = Inf, max_iter::Int = 10, improvement_tol::Float64 = 1.0e-3, termination::Symbol = :avg, relaxed::Bool = true, ext::Dict{Symbol,<:Any} = Dict{Symbol,Any}(:pump_breakpoints=>5), kwargs...) # Print a message with relevant algorithm limit information. @@ -195,6 +195,9 @@ function run_obbt_owf!(data::Dict{String,<:Any}, optimizer; # 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) From b808569ccb67ed1104cbc441e2c8ff2ed182623f Mon Sep 17 00:00:00 2001 From: Byron Tasseff Date: Tue, 18 Aug 2020 17:06:36 -0600 Subject: [PATCH 14/24] Initial version of reduced OBBT seems to work on the Van Zyl network. --- src/core/base.jl | 9 ++ src/core/constraint.jl | 6 +- src/core/constraint_template.jl | 10 +- src/core/ref.jl | 20 +++- src/core/variable.jl | 20 ++++ src/prob/wf.jl | 6 ++ src/util/obbt.jl | 164 ++++++++++++++++++++++++++++---- 7 files changed, 208 insertions(+), 27 deletions(-) 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 90e0acb6..7955a1f6 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,7 @@ 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 diff --git a/src/core/ref.jl b/src/core/ref.jl index 0d7893ae..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)) @@ -103,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(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) 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/prob/wf.jl b/src/prob/wf.jl index 5b9f6f47..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) @@ -63,8 +64,12 @@ function build_wf(wm::AbstractWaterModel) 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 end @@ -98,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) diff --git a/src/util/obbt.jl b/src/util/obbt.jl index 90319866..5a94fad9 100644 --- a/src/util/obbt.jl +++ b/src/util/obbt.jl @@ -68,7 +68,26 @@ function _get_edge_bound_dict(wm::AbstractWaterModel, nw::Int, comp::Symbol) end -function _create_modifications(wm::AbstractWaterModel) +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) for nw in sort(collect(nw_ids(wm))) @@ -109,7 +128,13 @@ 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]) - data_index = data["nw"][nw_str][comp_str][index_str] + + 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 @@ -120,52 +145,134 @@ function _update_modifications!(data::Dict{String,Any}, var_index_set::Array, bo 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 - 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] + + 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 _print_average_widths(data::Dict{String,<:Any}) +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(length(data["node"]))] 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 -function run_obbt_owf!(data::Dict{String,<:Any}, optimizer; +"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"]) + end + end + end + + for (i, tank) in ts_data["tank"] + tank["dispatchable"] = true + 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 +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, rel_gap_tol::Float64 = Inf, max_iter::Int = 10, - improvement_tol::Float64 = 1.0e-3, termination::Symbol = :avg, relaxed::Bool = true, + Inf, upper_bound_constraint::Bool = false, rel_gap_tol::Float64 = Inf, max_iter::Int = 100, + improvement_tol::Float64 = 1.0e-6, termination::Symbol = :avg, relaxed::Bool = true, 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, rel_gap_tol, upper_bound_constraint) # Instantiate the bound tightening model and relax integrality, if specified. - bt = instantiate_model(data, model_type, WaterModels.build_mn_owf; ext=ext) + bt = instantiate_model(data, model_type, build_type; ext=ext) upper_bound_constraint && _constraint_obj_bound(bt, upper_bound) - relaxed && _relax_model!(bt) + relaxed && _relax_model!(bt) # Relax integrality, if required. # Build the dictionary and sets that will store to the network. - modifications = _create_modifications(bt) + modifications = use_reduced_network ? _create_modifications_reduced(bt) : _create_modifications_mn(bt) var_index_set = vcat(_get_head_index_set(bt), _get_flow_index_set(bt)) bounds = [_get_existing_bounds(modifications, vid) for vid in var_index_set] - # Print the initial average bound widths. - _print_average_widths(modifications) + # 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 @@ -202,14 +309,31 @@ function run_obbt_owf!(data::Dict{String,<:Any}, optimizer; _update_modifications!(modifications, var_index_set, bounds) _IM.update_data!(data, modifications) - # Print the current average bound widths. - _print_average_widths(data) + # 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 - # Instantiate the new model and set the optimizer. - bt = instantiate_model(data, model_type, build_mn_owf; ext=ext) - upper_bound_constraint && _constraint_obj_bound(bt, upper_bound) - relaxed && _relax_model!(bt) - JuMP.set_optimizer(bt.model, optimizer) + if use_reduced_network + _revert_reduced_data!(data) end end From 9930580df02eb9a9cc8213280792f594deba9915 Mon Sep 17 00:00:00 2001 From: Byron Tasseff Date: Tue, 18 Aug 2020 21:08:26 -0600 Subject: [PATCH 15/24] Address reservoirs in network reduction. --- src/form/directed_flow.jl | 2 +- src/form/micp.jl | 4 ++-- src/form/milp.jl | 4 ++-- src/form/milpr.jl | 4 ++-- src/prob/owf.jl | 2 +- src/util/obbt.jl | 19 +++++++++++++++++-- 6 files changed, 25 insertions(+), 10 deletions(-) diff --git a/src/form/directed_flow.jl b/src/form/directed_flow.jl index d34d7f01..ae0da6bf 100644 --- a/src/form/directed_flow.jl +++ b/src/form/directed_flow.jl @@ -319,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) diff --git a/src/form/micp.jl b/src/form/micp.jl index 44256967..c8545681 100644 --- a/src/form/micp.jl +++ b/src/form/micp.jl @@ -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) @@ -153,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 395edbd7..254737b8 100644 --- a/src/form/milpr.jl +++ b/src/form/milpr.jl @@ -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) @@ -201,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/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/util/obbt.jl b/src/util/obbt.jl index 5a94fad9..081fb141 100644 --- a/src/util/obbt.jl +++ b/src/util/obbt.jl @@ -111,13 +111,13 @@ function _create_modifications_mn(wm::AbstractWaterModel) end -function _get_head_index_set(wm::AbstractWaterModel; width::Float64=0.1, prec::Int=4) +function _get_head_index_set(wm::AbstractWaterModel; width::Float64=1.0e-4, prec::Int=5) 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-5, prec::Int=7) +function _get_flow_index_set(wm::AbstractWaterModel; width::Float64=1.0e-6, prec::Int=7) 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]...) @@ -207,6 +207,9 @@ function _make_reduced_data!(ts_data::Dict{String,<:Any}) 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 @@ -214,6 +217,10 @@ function _make_reduced_data!(ts_data::Dict{String,<:Any}) for (i, tank) in ts_data["tank"] tank["dispatchable"] = true end + + for (i, reservoir) in ts_data["reservoir"] + reservoir["dispatchable"] = true + end end function _revert_reduced_data!(ts_data::Dict{String,<:Any}) @@ -233,6 +240,13 @@ function _revert_reduced_data!(ts_data::Dict{String,<:Any}) 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 @@ -265,6 +279,7 @@ function run_obbt_owf!(data::Dict{String,<:Any}, optimizer; use_reduced_network: # Build the dictionary and sets that will store to the network. modifications = use_reduced_network ? _create_modifications_reduced(bt) : _create_modifications_mn(bt) var_index_set = vcat(_get_head_index_set(bt), _get_flow_index_set(bt)) + #var_index_set = _get_flow_index_set(bt) bounds = [_get_existing_bounds(modifications, vid) for vid in var_index_set] # Get the initial average bound widths. From e22f7395a8c6a9369a37d931bdda5460f52b9896 Mon Sep 17 00:00:00 2001 From: Byron Tasseff Date: Tue, 18 Aug 2020 21:42:28 -0600 Subject: [PATCH 16/24] Modify widths and precision. --- src/util/obbt.jl | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/util/obbt.jl b/src/util/obbt.jl index 081fb141..0408ec3d 100644 --- a/src/util/obbt.jl +++ b/src/util/obbt.jl @@ -87,6 +87,7 @@ function _create_modifications_reduced(wm::AbstractWaterModel) return data end + function _create_modifications_mn(wm::AbstractWaterModel) data = Dict{String,Any}("nw"=>Dict{String,Any}(), "per_unit"=>false) @@ -111,13 +112,13 @@ function _create_modifications_mn(wm::AbstractWaterModel) end -function _get_head_index_set(wm::AbstractWaterModel; width::Float64=1.0e-4, prec::Int=5) +function _get_head_index_set(wm::AbstractWaterModel; width::Float64=1.0e-6, prec::Int=10) 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-6, prec::Int=7) +function _get_flow_index_set(wm::AbstractWaterModel; width::Float64=1.0e-8, prec::Int=12) 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]...) @@ -223,6 +224,7 @@ function _make_reduced_data!(ts_data::Dict{String,<:Any}) 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. From bf87e9fa078cf14e9451282add19b9287f36cda2 Mon Sep 17 00:00:00 2001 From: Byron Tasseff Date: Tue, 18 Aug 2020 21:55:36 -0600 Subject: [PATCH 17/24] Fix dispatchable reservoirs. --- src/util/obbt.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/util/obbt.jl b/src/util/obbt.jl index 0408ec3d..7c58db44 100644 --- a/src/util/obbt.jl +++ b/src/util/obbt.jl @@ -220,7 +220,9 @@ function _make_reduced_data!(ts_data::Dict{String,<:Any}) end for (i, reservoir) in ts_data["reservoir"] - reservoir["dispatchable"] = true + if string(reservoir["node"]) in keys(ts_data["time_series"]["node"]) + reservoir["dispatchable"] = true + end end end @@ -281,7 +283,6 @@ function run_obbt_owf!(data::Dict{String,<:Any}, optimizer; use_reduced_network: # Build the dictionary and sets that will store to the network. modifications = use_reduced_network ? _create_modifications_reduced(bt) : _create_modifications_mn(bt) var_index_set = vcat(_get_head_index_set(bt), _get_flow_index_set(bt)) - #var_index_set = _get_flow_index_set(bt) bounds = [_get_existing_bounds(modifications, vid) for vid in var_index_set] # Get the initial average bound widths. From b3438db8ce9d551a2b8d18cb83bd3dbb412d5124 Mon Sep 17 00:00:00 2001 From: Byron Tasseff Date: Tue, 18 Aug 2020 21:59:14 -0600 Subject: [PATCH 18/24] Fix dispatchable reservoirs. --- src/util/obbt.jl | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/util/obbt.jl b/src/util/obbt.jl index 7c58db44..b2d443d5 100644 --- a/src/util/obbt.jl +++ b/src/util/obbt.jl @@ -220,8 +220,10 @@ function _make_reduced_data!(ts_data::Dict{String,<:Any}) end for (i, reservoir) in ts_data["reservoir"] - if string(reservoir["node"]) in keys(ts_data["time_series"]["node"]) - reservoir["dispatchable"] = true + 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 From 2e697a9b687fa8581ca377d04b2cd4004955e76c Mon Sep 17 00:00:00 2001 From: Byron Tasseff Date: Thu, 20 Aug 2020 10:34:50 -0600 Subject: [PATCH 19/24] Add a function for making tank volumes dispatchable at the first time step. --- src/core/data.jl | 12 ++++++++++++ 1 file changed, 12 insertions(+) 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)]) From 7a49c3e8047163f67d0a0a42fa62359d38932941 Mon Sep 17 00:00:00 2001 From: Byron Tasseff Date: Thu, 20 Aug 2020 10:49:34 -0600 Subject: [PATCH 20/24] Correct integration after a dispatchable tank state. --- src/core/constraint_template.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/core/constraint_template.jl b/src/core/constraint_template.jl index 7955a1f6..3e44f060 100644 --- a/src/core/constraint_template.jl +++ b/src/core/constraint_template.jl @@ -189,7 +189,7 @@ function constraint_tank_state(wm::AbstractWaterModel, i::Int, nw_1::Int, nw_2:: end # Only set the tank state if the tank is nondispatchable. - if !ref(wm, nw_1, :tank, i)["dispatchable"] + 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) From fc9015bc885c01361b456ecd5f570a3c9ee5be2d Mon Sep 17 00:00:00 2001 From: Byron Tasseff Date: Thu, 20 Aug 2020 16:08:07 -0600 Subject: [PATCH 21/24] Change default precisions and bound widths. --- src/util/obbt.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/util/obbt.jl b/src/util/obbt.jl index b2d443d5..245bc159 100644 --- a/src/util/obbt.jl +++ b/src/util/obbt.jl @@ -112,13 +112,13 @@ function _create_modifications_mn(wm::AbstractWaterModel) end -function _get_head_index_set(wm::AbstractWaterModel; width::Float64=1.0e-6, prec::Int=10) +function _get_head_index_set(wm::AbstractWaterModel; width::Float64=1.0e-2, prec::Int=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-8, prec::Int=12) +function _get_flow_index_set(wm::AbstractWaterModel; width::Float64=1.0e-5, prec::Int=7) 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]...) From 7835dd45a7dc84e3606667538507f1de0c6cb946 Mon Sep 17 00:00:00 2001 From: Byron Tasseff Date: Fri, 21 Aug 2020 09:56:37 -0600 Subject: [PATCH 22/24] Add OBBT utility tests. --- src/util/obbt.jl | 6 +++--- test/runtests.jl | 2 ++ test/util.jl | 14 ++++++++++++++ 3 files changed, 19 insertions(+), 3 deletions(-) create mode 100644 test/util.jl diff --git a/src/util/obbt.jl b/src/util/obbt.jl index 245bc159..33b67299 100644 --- a/src/util/obbt.jl +++ b/src/util/obbt.jl @@ -1,5 +1,5 @@ ### Optimization-based Bound Tightening for Optimal Water Flow Problems ### -### This is built upon https://github.com/lanl-ansi/WaterModels.jl/blob/master/src/util/obbt.jl +### 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. @@ -89,7 +89,7 @@ end function _create_modifications_mn(wm::AbstractWaterModel) - data = Dict{String,Any}("nw"=>Dict{String,Any}(), "per_unit"=>false) + 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}() @@ -181,7 +181,7 @@ 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(length(data["node"]))] + 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 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 From b15aab04acac4b1ff260140823186dfcb30409d6 Mon Sep 17 00:00:00 2001 From: Byron Tasseff Date: Fri, 21 Aug 2020 10:35:31 -0600 Subject: [PATCH 23/24] Update CHANGELOG. --- CHANGELOG.md | 1 + 1 file changed, 1 insertion(+) 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 From 86d796bfb78047395bb69c5717d3b6c1f28267f5 Mon Sep 17 00:00:00 2001 From: Byron Tasseff Date: Fri, 21 Aug 2020 11:01:47 -0600 Subject: [PATCH 24/24] Address @ccoffrin's comments, tidy up docstrings. --- .travis.yml | 4 ++-- src/util/obbt.jl | 46 +++++++++++++++++----------------------------- 2 files changed, 19 insertions(+), 31 deletions(-) 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/src/util/obbt.jl b/src/util/obbt.jl index 33b67299..b9033047 100644 --- a/src/util/obbt.jl +++ b/src/util/obbt.jl @@ -6,8 +6,6 @@ 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. -Interested readers are referred to the paper "Strengthening Convex Relaxations with Bound -Tightening for Water Network Optimization." # Example @@ -15,29 +13,21 @@ The function can be invoked as follows: ``` ipopt = JuMP.optimizer_with_attributes(Ipopt.Optimizer) -data, stats = run_obbt_owf!("examples/data/epanet/van_zyl.inp", ipopt) +run_obbt_owf!("examples/data/epanet/van_zyl.inp", ipopt) ``` -`data` contains the parsed network data with tightened bounds. -`stats` contains information output from the bound tightening algorithm. - # 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. -* `rel_gap_tol`: tolerance used to terminate the algorithm when the objective - value of the relaxation is close to the upper bound specified using the - `upper_bound` keyword. -* `min_width`: domain beyond which bound tightening is not performed. -* `termination`: Bound-tightening algorithm terminates if the improvement in - the average or maximum bound improvement, specified using either the - `termination = :avg` or the `termination = :max` option, is less than - `improvement_tol`. -* `precision`: number of decimal digits to round the tightened bounds to. +* `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) @@ -112,13 +102,13 @@ function _create_modifications_mn(wm::AbstractWaterModel) end -function _get_head_index_set(wm::AbstractWaterModel; width::Float64=1.0e-2, prec::Int=4) +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-5, prec::Int=7) +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]...) @@ -258,8 +248,8 @@ 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, rel_gap_tol::Float64 = Inf, max_iter::Int = 100, - improvement_tol::Float64 = 1.0e-6, termination::Symbol = :avg, relaxed::Bool = true, + 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.") @@ -275,7 +265,7 @@ function run_obbt_owf!(data::Dict{String,<:Any}, optimizer; use_reduced_network: end # Check for keyword argument inconsistencies. - _check_obbt_options(upper_bound, rel_gap_tol, upper_bound_constraint) + _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) @@ -284,7 +274,9 @@ function run_obbt_owf!(data::Dict{String,<:Any}, optimizer; use_reduced_network: # Build the dictionary and sets that will store to the network. modifications = use_reduced_network ? _create_modifications_reduced(bt) : _create_modifications_mn(bt) - var_index_set = vcat(_get_head_index_set(bt), _get_flow_index_set(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. @@ -358,21 +350,17 @@ function run_obbt_owf!(data::Dict{String,<:Any}, optimizer; use_reduced_network: end -function _check_obbt_options(ub::Float64, rel_gap::Float64, ub_constraint::Bool) +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 - - if !isinf(rel_gap) && isinf(ub) - Memento.error(_LOGGER, "[OBBT] The option \"rel_gap_tol\" is specified without providing an upper bound.") - end end -function _compute_corrected_bounds(lb::Float64, ub::Float64, lb_old::Float64, ub_old::Float64, width::Float64, prec::Int) +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(10.0^prec * lb) * inv(10.0^prec), lb_old) - ub = min(ceil(10.0^prec * ub) * inv(10.0^prec), ub_old) + 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