From 8cc94dc1abaf66c65cc53a9855944e426de894c9 Mon Sep 17 00:00:00 2001 From: Byron Tasseff Date: Sat, 22 Aug 2020 12:19:28 -0600 Subject: [PATCH 1/3] Explicitly separate fixed pipes from design pipes. --- docs/make.jl | 2 +- docs/src/index.md | 2 +- docs/src/quickguide.md | 22 +-- {test => examples}/data/json/shamir.json | 0 src/core/base.jl | 8 +- src/core/constraint.jl | 20 +- src/core/constraint_template.jl | 45 +++-- src/core/data.jl | 4 +- src/core/objective.jl | 2 +- src/core/ref.jl | 2 +- src/core/variable.jl | 6 +- src/form/directed_flow.jl | 224 +++++++++++------------ src/form/micp.jl | 4 +- src/form/milp.jl | 42 ++--- src/form/milpr.jl | 4 +- src/form/nlp.jl | 4 +- src/form/undirected_flow.jl | 75 ++++---- src/prob/des.jl | 36 ++-- test/des.jl | 2 +- test/io.jl | 10 +- test/warm_start.jl | 2 +- 21 files changed, 266 insertions(+), 250 deletions(-) rename {test => examples}/data/json/shamir.json (100%) diff --git a/docs/make.jl b/docs/make.jl index 5382dc9d..8c87cec7 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -2,7 +2,7 @@ using Documenter, WaterModels makedocs( modules = [WaterModels], - format = Documenter.HTML(analytics="UA-367975-10", mathengine=Documenter.MathJax()), + format = Documenter.HTML(analytics="UA-367975-10", mathengine=Documenter.MathJax(), prettyurls=false), sitename = "WaterModels", authors = "Byron Tasseff and contributors", pages = [ diff --git a/docs/src/index.md b/docs/src/index.md index 49c786e0..ed013329 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -7,7 +7,7 @@ CurrentModule = WaterModels ## Overview WaterModels.jl is a Julia/JuMP package for steady state potable water distribution network optimization. It is designed to enable the computational evaluation of historical and emerging water network optimization formulations and algorithms using a common platform. -The code is engineered to decouple [Problem Specifications](@ref) (e.g., water flow, optimal water flow, network design) from [Network Formulations](@ref) (e.g., mixed-integer linear, mixed-integer nonlinear). +The code is specifically engineered to decouple [Problem Specifications](@ref) (e.g., water flow, optimal water flow, network design) from [Network Formulations](@ref) (e.g., mixed-integer linear, mixed-integer nonlinear). This decoupling enables the definition of a wide variety of water network optimization formulations and their comparison across several common problem specifications. ## Installation diff --git a/docs/src/quickguide.md b/docs/src/quickguide.md index 679d29a0..f9d07dbc 100644 --- a/docs/src/quickguide.md +++ b/docs/src/quickguide.md @@ -15,18 +15,13 @@ For the current development version, install the package using ] add WaterModels#master ``` -Next, test that the package works by executing +Finally, test that the package works as expected by executing ```julia ] test WaterModels ``` -Install [Cbc](https://github.com/JuliaOpt/Cbc.jl) using -```julia -] add Cbc -``` - ## Solving a Network Design Problem -Once the above dependencies have been installed, obtain the files [`shamir.inp`](https://raw.githubusercontent.com/lanl-ansi/WaterModels.jl/master/test/data/epanet/shamir.inp) and [`shamir.json`](https://raw.githubusercontent.com/lanl-ansi/WaterModels.jl/master/test/data/json/shamir.json). +Once the above dependencies have been installed, obtain the files [`shamir.inp`](https://raw.githubusercontent.com/lanl-ansi/WaterModels.jl/master/examples/data/epanet/shamir.inp) and [`shamir.json`](https://raw.githubusercontent.com/lanl-ansi/WaterModels.jl/master/examples/data/json/shamir.json). Here, `shamir.inp` is an EPANET file describing a simple seven-node, eight-link water distribution network with one reservoir, six junctions, and eight pipes. In accord, `shamir.json` is a JSON file specifying possible pipe diameters and associated costs per unit length, per diameter setting. The combination of data from these two files provides the required information to set up a corresponding network design problem, where the goal is to select the most cost-efficient pipe diameters while satisfying all demand in the network. @@ -36,12 +31,12 @@ To read in the EPANET data, execute the following: ```julia using WaterModels -data = parse_file("shamir.inp") +data = parse_file("examples/data/epanet/shamir.inp") ``` Next, to read in the JSON diameter and cost data, execute the following: ```julia -modifications = parse_file("shamir.json") +modifications = parse_file("examples/data/json/shamir.json") ``` To merge these data together, InfrastructureModels is used to update the original network data using @@ -52,7 +47,6 @@ WaterModels._IM.update_data!(data, modifications) Finally, the MILP formulation for the network design specification can be solved using ```julia import Cbc - run_des(data, MILPWaterModel, Cbc.Optimizer) ``` @@ -67,7 +61,7 @@ However, because of the finer discretization, a better approximation of the phys Instead of linear approximation, head loss curves can also be linearly outer-approximated via the MILP-R formulation. This formulation employs less strict requirements and avoids the use of binary variables, but solutions (e.g., diameters) may not necessarily be feasible with respect to the full (nonconvex) water network physics. -To employ five outer approximation points per (positive or negative) head loss curve in this formulation, the following can be executed +To employ five outer approximation points per (positive or negative) head loss curve in this formulation, the following can be executed: ```julia run_des(data, MILPRWaterModel, Cbc.Optimizer, ext=Dict(:pipe_breakpoints=>5)) ``` @@ -88,7 +82,7 @@ result["objective"] The `"solution"` field contains detailed information about the solution produced by the `run` method. For example, the following dictionary comprehension can be used to inspect the flows in the solution: ``` -Dict(name => data["q"] for (name, data) in result["solution"]["link"]) +Dict(name => data["q"] for (name, data) in result["solution"]["pipe_des"]) ``` For more information about WaterModels result data see the [WaterModels Result Data Format](@ref) section. @@ -99,11 +93,9 @@ Mixed-integer nonconvex formulations can be solved with dedicated solvers, as we For example, the full mixed-integer nonconvex formulation for design (NLP) can be solved via ```julia import KNITRO - -run_des(data, NCNLPWaterModel, KNITRO.Optimizer) +run_des(data, NLPWaterModel, KNITRO.Optimizer) ``` and the mixed-integer convex formulation (MICP) can be solved via - ```julia run_des(data, MICPWaterModel, KNITRO.Optimizer) ``` diff --git a/test/data/json/shamir.json b/examples/data/json/shamir.json similarity index 100% rename from test/data/json/shamir.json rename to examples/data/json/shamir.json diff --git a/src/core/base.jl b/src/core/base.jl index b39e555c..f9a73988 100644 --- a/src/core/base.jl +++ b/src/core/base.jl @@ -157,10 +157,6 @@ function _ref_add_core!(nw_refs::Dict{Int,<:Any}) # 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]) @@ -172,8 +168,8 @@ function _ref_add_core!(nw_refs::Dict{Int,<:Any}) ref[:shutoff_valve] = filter(x -> x.second["has_shutoff_valve"], ref[:pipe]) ref[:pipe] = filter(x -> !x.second["has_check_valve"], ref[:pipe]) ref[:pipe] = filter(x -> !x.second["has_shutoff_valve"], ref[:pipe]) - ref[:pipe_des] = filter(is_des_pipe, ref[:pipe]) - ref[:pipe_fixed] = filter(!is_des_pipe, ref[:pipe]) + ref[:des_pipe] = filter(is_des_pipe, ref[:pipe]) + ref[:pipe] = filter(!is_des_pipe, ref[:pipe]) # Create mappings of "from" and "to" arcs for link- (i.e., edge-) type components. for name in diff --git a/src/core/constraint.jl b/src/core/constraint.jl index 27cef9b4..a14451e9 100644 --- a/src/core/constraint.jl +++ b/src/core/constraint.jl @@ -15,21 +15,24 @@ end """ constraint_flow_conservation( - wm, n, 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, demand) + wm, n, i, check_valve_fr, check_valve_to, pipe_fr, pipe_to, des_pipe_fr, + des_pipe_to, pump_fr, pump_to, pressure_reducing_valve_fr, + pressure_reducing_valve_to, shutoff_valve_fr, shutoff_valve_to, reservoirs, tanks, + dispatachable_junctions, fixed_demand) """ function constraint_flow_conservation( wm::AbstractWaterModel, n::Int, i::Int, check_valve_fr::Array{Int64,1}, check_valve_to::Array{Int64,1}, pipe_fr::Array{Int64,1}, pipe_to::Array{Int64,1}, - 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}, + des_pipe_fr::Array{Int64,1}, des_pipe_to::Array{Int64,1}, pump_fr::Array{Int64,1}, + pump_to::Array{Int64,1}, pressure_reducing_valve_fr::Array{Int64,1}, + pressure_reducing_valve_to::Array{Int64,1}, shutoff_valve_fr::Array{Int64,1}, + shutoff_valve_to::Array{Int64,1}, reservoirs::Array{Int64,1}, tanks::Array{Int64,1}, dispatchable_junctions::Array{Int64,1}, fixed_demand::Float64) # 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) + q_pipe = var(wm, n, :q_pipe) + q_des_pipe = var(wm, n, :q_des_pipe_sum) + q_pump = var(wm, n, :q_pump) q_pressure_reducing_valve = var(wm, n, :q_pressure_reducing_valve) q_shutoff_valve = var(wm, n, :q_shutoff_valve) @@ -38,6 +41,7 @@ function constraint_flow_conservation( sum(q_check_valve[a] for a in check_valve_fr) + sum(q_check_valve[a] for a in check_valve_to) - sum(q_pipe[a] for a in pipe_fr) + sum(q_pipe[a] for a in pipe_to) - + sum(q_des_pipe[a] for a in des_pipe_fr) + sum(q_des_pipe[a] for a in des_pipe_to) - sum(q_pump[a] for a in pump_fr) + sum(q_pump[a] for a in pump_to) - sum(q_pressure_reducing_valve[a] for a in pressure_reducing_valve_fr) + sum(q_pressure_reducing_valve[a] for a in pressure_reducing_valve_to) - diff --git a/src/core/constraint_template.jl b/src/core/constraint_template.jl index 3e44f060..6b7513cb 100644 --- a/src/core/constraint_template.jl +++ b/src/core/constraint_template.jl @@ -32,6 +32,8 @@ function constraint_flow_conservation(wm::AbstractWaterModel, i::Int; nw::Int=wm check_valve_to = _collect_comps_to(wm, i, :check_valve; nw=nw) pipe_fr = _collect_comps_fr(wm, i, :pipe; nw=nw) pipe_to = _collect_comps_to(wm, i, :pipe; nw=nw) + des_pipe_fr = _collect_comps_fr(wm, i, :des_pipe; nw=nw) + des_pipe_to = _collect_comps_to(wm, i, :des_pipe; nw=nw) pump_fr = _collect_comps_fr(wm, i, :pump; nw=nw) pump_to = _collect_comps_to(wm, i, :pump; nw=nw) pressure_reducing_valve_fr = _collect_comps_fr(wm, i, :pressure_reducing_valve; nw=nw) @@ -57,9 +59,10 @@ function constraint_flow_conservation(wm::AbstractWaterModel, i::Int; nw::Int=wm # Add the flow conservation constraint. 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, dispatchable_junctions, net_fixed_demand) + wm, nw, i, check_valve_fr, check_valve_to, pipe_fr, pipe_to, des_pipe_fr, + des_pipe_to, pump_fr, pump_to, pressure_reducing_valve_fr, + pressure_reducing_valve_to, shutoff_valve_fr, shutoff_valve_to, reservoirs, tanks, + dispatchable_junctions, net_fixed_demand) end @@ -69,6 +72,8 @@ function constraint_node_directionality(wm::AbstractWaterModel, i::Int; nw::Int= check_valve_to = _collect_comps_to(wm, i, :check_valve; nw=nw) pipe_fr = _collect_comps_fr(wm, i, :pipe; nw=nw) pipe_to = _collect_comps_to(wm, i, :pipe; nw=nw) + des_pipe_fr = _collect_comps_fr(wm, i, :des_pipe; nw=nw) + des_pipe_to = _collect_comps_to(wm, i, :des_pipe; nw=nw) pump_fr = _collect_comps_fr(wm, i, :pump; nw=nw) pump_to = _collect_comps_to(wm, i, :pump; nw=nw) pressure_reducing_valve_fr = _collect_comps_fr(wm, i, :pressure_reducing_valve; nw=nw) @@ -83,12 +88,12 @@ function constraint_node_directionality(wm::AbstractWaterModel, i::Int; nw::Int= num_components = length(junctions) + length(tanks) + length(reservoirs) # Get the in degree of node `i`. - in_length = length(check_valve_to) + length(pipe_to) + length(pump_to) + - length(pressure_reducing_valve_to) + length(shutoff_valve_to) + in_length = length(check_valve_to) + length(pipe_to) + length(des_pipe_to) + + length(pump_to) + length(pressure_reducing_valve_to) + length(shutoff_valve_to) # 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) + out_length = length(check_valve_fr) + length(pipe_fr) + length(des_pipe_fr) + + length(pump_fr) + length(pressure_reducing_valve_fr) + length(shutoff_valve_fr) # Check if node directionality constraints should be added. if num_components == 0 && in_length + out_length == 2 @@ -97,9 +102,9 @@ function constraint_node_directionality(wm::AbstractWaterModel, i::Int; nw::Int= # 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) + wm, nw, i, check_valve_fr, check_valve_to, pipe_fr, pipe_to, des_pipe_fr, + des_pipe_to, pump_fr, pump_to, pressure_reducing_valve_fr, + pressure_reducing_valve_to, shutoff_valve_fr, shutoff_valve_to) end end @@ -111,6 +116,8 @@ function constraint_sink_directionality(wm::AbstractWaterModel, i::Int; nw::Int= check_valve_to = _collect_comps_to(wm, i, :check_valve; nw=nw) pipe_fr = _collect_comps_fr(wm, i, :pipe; nw=nw) pipe_to = _collect_comps_to(wm, i, :pipe; nw=nw) + des_pipe_fr = _collect_comps_fr(wm, i, :des_pipe; nw=nw) + des_pipe_to = _collect_comps_to(wm, i, :des_pipe; nw=nw) pump_fr = _collect_comps_fr(wm, i, :pump; nw=nw) pump_to = _collect_comps_to(wm, i, :pump; nw=nw) pressure_reducing_valve_fr = _collect_comps_fr(wm, i, :pressure_reducing_valve; nw=nw) @@ -123,9 +130,9 @@ function constraint_sink_directionality(wm::AbstractWaterModel, i::Int; nw::Int= # Add the sink flow directionality constraint. constraint_sink_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) + wm, nw, i, check_valve_fr, check_valve_to, pipe_fr, pipe_to, des_pipe_fr, + des_pipe_to, pump_fr, pump_to, pressure_reducing_valve_fr, + pressure_reducing_valve_to, shutoff_valve_fr, shutoff_valve_to) end @@ -147,6 +154,8 @@ function constraint_source_directionality(wm::AbstractWaterModel, i::Int; nw::In check_valve_to = _collect_comps_to(wm, i, :check_valve; nw=nw) pipe_fr = _collect_comps_fr(wm, i, :pipe; nw=nw) pipe_to = _collect_comps_to(wm, i, :pipe; nw=nw) + des_pipe_fr = _collect_comps_fr(wm, i, :des_pipe; nw=nw) + des_pipe_to = _collect_comps_to(wm, i, :des_pipe; nw=nw) pump_fr = _collect_comps_fr(wm, i, :pump; nw=nw) pump_to = _collect_comps_to(wm, i, :pump; nw=nw) pressure_reducing_valve_fr = _collect_comps_fr(wm, i, :pressure_reducing_valve; nw=nw) @@ -159,9 +168,9 @@ function constraint_source_directionality(wm::AbstractWaterModel, i::Int; nw::In # Add the source flow directionality constraint. constraint_source_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) + wm, nw, i, check_valve_fr, check_valve_to, pipe_fr, pipe_to, des_pipe_fr, + des_pipe_to, pump_fr, pump_to, pressure_reducing_valve_fr, + pressure_reducing_valve_to, shutoff_valve_fr, shutoff_valve_to) end @@ -219,7 +228,7 @@ end function constraint_pipe_head_loss_des(wm::AbstractWaterModel, a::Int; nw::Int=wm.cnw, kwargs...) alpha = ref(wm, nw, :alpha) - pipe = ref(wm, nw, :pipe, a) + pipe = ref(wm, nw, :des_pipe, a) resistances = ref(wm, nw, :resistance, a) i, j = pipe["node_fr"], pipe["node_to"] @@ -246,7 +255,7 @@ end function constraint_pipe_head_loss_ub_des(wm::AbstractWaterModel, a::Int; nw::Int=wm.cnw) - alpha, L = [ref(wm, nw, :alpha), ref(wm, nw, :pipe, a)["length"]] + alpha, L = ref(wm, nw, :alpha), ref(wm, nw, :des_pipe, a)["length"] pipe_resistances = ref(wm, nw, :resistance, a) constraint_pipe_head_loss_ub_des(wm, nw, a, alpha, L, pipe_resistances) end diff --git a/src/core/data.jl b/src/core/data.jl index a0b977fc..03ee25a8 100644 --- a/src/core/data.jl +++ b/src/core/data.jl @@ -257,9 +257,9 @@ function set_start_undirected_flow_rate_des!(data::Dict{String, <:Any}) for (a, pipe) in filter(is_des_pipe, data["pipe"]) num_resistances = length(resistances[a]) - pipe["q_des_start"] = zeros(Float64, num_resistances) + pipe["q_des_pipe_start"] = zeros(Float64, num_resistances) r_id, val = findmax(pipe["x_res_start"]) - pipe["q_des_start"][r_id] = pipe["q"] + pipe["q_des_pipe_start"][r_id] = pipe["q"] end end diff --git a/src/core/objective.jl b/src/core/objective.jl index d857c55f..28106d02 100644 --- a/src/core/objective.jl +++ b/src/core/objective.jl @@ -23,7 +23,7 @@ function objective_des(wm::AbstractWaterModel) obj = JuMP.AffExpr(0.0) for n in nw_ids(wm) - for (a, pipe) in ref(wm, n, :pipe_des) + for (a, pipe) in ref(wm, n, :des_pipe) costs = pipe["length"] .* ref(wm, n, :resistance_cost, a) JuMP.add_to_expression!(obj, sum(costs .* var(wm, n, :x_res, a))) end diff --git a/src/core/ref.jl b/src/core/ref.jl index ef53ae8c..3fed696c 100644 --- a/src/core/ref.jl +++ b/src/core/ref.jl @@ -130,7 +130,7 @@ function calc_flow_bounds(wm::AbstractWaterModel, n::Int=wm.cnw) lb, ub = Dict{String,Any}(), Dict{String,Any}() - for name in ["check_valve", "pipe", "shutoff_valve"] + for name in ["check_valve", "pipe", "des_pipe", "shutoff_valve"] lb[name], ub[name] = Dict{Int,Array{Float64}}(), Dict{Int,Array{Float64}}() for (a, comp) in ref(wm, n, Symbol(name)) diff --git a/src/core/variable.jl b/src/core/variable.jl index 24223da7..47737c2d 100644 --- a/src/core/variable.jl +++ b/src/core/variable.jl @@ -199,12 +199,12 @@ where one denotes that the given resistance is active in the design." function variable_resistance(wm::AbstractWaterModel; nw::Int=wm.cnw, report::Bool=true) x_res = var(wm, nw)[:x_res] = Dict{Int, Array{JuMP.VariableRef}}() - for a in ids(wm, nw, :pipe_des) + for a in ids(wm, nw, :des_pipe) n_r = length(ref(wm, nw, :resistance, a)) # Number of resistances. var(wm, nw, :x_res)[a] = JuMP.@variable(wm.model, [r in 1:n_r], binary=true, base_name="$(nw)_x_res[$(a)]", - start=comp_start_value(ref(wm, nw, :pipe_des, a), "x_res_start", r)) + start=comp_start_value(ref(wm, nw, :des_pipe, a), "x_res_start", r)) end - report && sol_component_value(wm, nw, :pipe_des, :x_res, ids(wm, nw, :pipe_des), x_res) + report && sol_component_value(wm, nw, :des_pipe, :x_res, ids(wm, nw, :des_pipe), x_res) end diff --git a/src/form/directed_flow.jl b/src/form/directed_flow.jl index ae0da6bf..a5222a2b 100644 --- a/src/form/directed_flow.jl +++ b/src/form/directed_flow.jl @@ -4,13 +4,6 @@ # qn is nonzero, qp should be zero. -function _collect_directions(wm::AbstractWaterModel, n::Int, a::Int) - comps = ["check_valve", "pipe", "pressure_reducing_valve", "pump", "shutoff_valve"] - comps = filter(x -> a in ids(wm, n, Symbol(x)), comps) - return sum(var(wm, n, Symbol("y_" * c), a) for c in comps) -end - - "Initialize variables associated with flow direction. If this variable is equal to one, the flow direction is from i to j. If it is equal to zero, the flow direction is from j to i." @@ -119,46 +112,55 @@ function variable_flow(wm::AbstractDirectedModel; nw::Int=wm.cnw, bounded::Bool= # Create directed flow binary direction variables (`y`) for each component. _variable_component_direction(wm, name; nw=nw, report=report) end + + # Create flow-related variables for design components. + variable_flow_des(wm; nw=nw, bounded=bounded, report=report) end "Create network design flow variables for directed flow formulations." function variable_flow_des(wm::AbstractDirectedModel; nw::Int=wm.cnw, bounded::Bool=true, report::Bool=true) - # Create dictionary for undirected design flow variables (qp_des and qn_des). - qp_des = var(wm, nw)[:qp_des] = Dict{Int,Array{JuMP.VariableRef}}() - qn_des = var(wm, nw)[:qn_des] = Dict{Int,Array{JuMP.VariableRef}}() + # Create dictionary for undirected design flow variables (qp_des_pipe and qn_des_pipe). + qp_des_pipe = var(wm, nw)[:qp_des_pipe] = Dict{Int,Array{JuMP.VariableRef}}() + qn_des_pipe = var(wm, nw)[:qn_des_pipe] = Dict{Int,Array{JuMP.VariableRef}}() # Initialize the variables. (The default start value of _q_eps is crucial.) - for a in ids(wm, nw, :pipe_des) - var(wm, nw, :qp_des)[a] = JuMP.@variable(wm.model, + for a in ids(wm, nw, :des_pipe) + var(wm, nw, :qp_des_pipe)[a] = JuMP.@variable(wm.model, [r in 1:length(ref(wm, nw, :resistance, a))], lower_bound=0.0, - base_name="$(nw)_qp_des[$(a)]", - start=comp_start_value(ref(wm, nw, :pipe_des, a), "qp_des_start", r, _q_eps)) + base_name="$(nw)_qp_des_pipe[$(a)]", + start=comp_start_value(ref(wm, nw, :des_pipe, a), "qp_des_pipe_start", r, _q_eps)) - var(wm, nw, :qn_des)[a] = JuMP.@variable(wm.model, + var(wm, nw, :qn_des_pipe)[a] = JuMP.@variable(wm.model, [r in 1:length(ref(wm, nw, :resistance, a))], lower_bound=0.0, - base_name="$(nw)_qn_des[$(a)]", - start=comp_start_value(ref(wm, nw, :pipe_des, a), "qn_des_start", r, _q_eps)) + base_name="$(nw)_qn_des_pipe[$(a)]", + start=comp_start_value(ref(wm, nw, :des_pipe, a), "qn_des_pipe_start", r, _q_eps)) end if bounded # If the variables are bounded, apply the bounds. q_lb, q_ub = calc_flow_bounds(wm, nw) - for a in ids(wm, nw, :pipe_des) + for a in ids(wm, nw, :des_pipe) for r in 1:length(ref(wm, nw, :resistance, a)) - JuMP.set_upper_bound(qp_des[a][r], max(0.0, q_ub["pipe"][a][r])) - JuMP.set_upper_bound(qn_des[a][r], max(0.0, -q_lb["pipe"][a][r])) + JuMP.set_upper_bound(qp_des_pipe[a][r], max(0.0, q_ub["des_pipe"][a][r])) + JuMP.set_upper_bound(qn_des_pipe[a][r], max(0.0, -q_lb["des_pipe"][a][r])) end end end - # Create expressions capturing the relationships among q, qp_des, and qn_des. - q = var(wm, nw)[:q] = JuMP.@expression( - wm.model, [a in ids(wm, nw, :pipe_des)], - sum(var(wm, nw, :qp_des, a)) - sum(var(wm, nw, :qn_des, a))) + # Create directed head difference (`dhp` and `dhn`) variables for each component. + _variable_component_head_difference(wm, "des_pipe"; nw=nw, bounded=bounded, report=report) + + # Create directed flow binary direction variables (`y`) for each component. + _variable_component_direction(wm, "des_pipe"; nw=nw, report=report) + + # Create expressions capturing the relationships among q, qp_des_pipe, and qn_des_pipe. + q = var(wm, nw)[:q_des_pipe_sum] = JuMP.@expression( + wm.model, [a in ids(wm, nw, :des_pipe)], + sum(var(wm, nw, :qp_des_pipe, a)) - sum(var(wm, nw, :qn_des_pipe, a))) # Initialize the solution reporting data structures. - report && sol_component_value(wm, nw, :pipe_des, :q, ids(wm, nw, :pipe_des), q) + report && sol_component_value(wm, nw, :des_pipe, :q, ids(wm, nw, :des_pipe), q) # Create resistance binary variables. variable_resistance(wm, nw=nw) @@ -362,10 +364,10 @@ end function constraint_flow_direction_selection_des(wm::AbstractDirectedModel, n::Int, a::Int, pipe_resistances) - y = var(wm, n, :y_pipe, a) + y = var(wm, n, :y_des_pipe, a) for r_id in 1:length(pipe_resistances) - qp, qn = var(wm, n, :qp_des, a)[r_id], var(wm, n, :qn_des, a)[r_id] + qp, qn = var(wm, n, :qp_des_pipe, a)[r_id], var(wm, n, :qn_des_pipe, a)[r_id] c_p = JuMP.@constraint(wm.model, qp <= JuMP.upper_bound(qp) * y) c_n = JuMP.@constraint(wm.model, qn <= JuMP.upper_bound(qn) * (1.0 - y)) append!(con(wm, n, :head_loss)[a], [c_p, c_n]) @@ -410,17 +412,17 @@ end function constraint_pipe_head_loss_ub_des(wm::AbstractDirectedModel, n::Int, a::Int, alpha::Float64, L::Float64, pipe_resistances) - dhp = var(wm, n, :dhp_pipe, a) - qp_des = var(wm, n, :qp_des, a) - qp_des_ub = JuMP.upper_bound.(qp_des) - slopes_p = pipe_resistances .* qp_des_ub.^(alpha - 1.0) - c_p = JuMP.@constraint(wm.model, inv(L)*dhp <= sum(slopes_p .* qp_des)) - - dhn = var(wm, n, :dhn_pipe, a) - qn_des = var(wm, n, :qn_des, a) - qn_des_ub = JuMP.upper_bound.(qn_des) - slopes_n = pipe_resistances .* qn_des_ub.^(alpha - 1.0) - c_n = JuMP.@constraint(wm.model, inv(L)*dhn <= sum(slopes_n .* qn_des)) + dhp = var(wm, n, :dhp_des_pipe, a) + qp_des_pipe = var(wm, n, :qp_des_pipe, a) + qp_des_pipe_ub = JuMP.upper_bound.(qp_des_pipe) + slopes_p = pipe_resistances .* qp_des_pipe_ub.^(alpha - 1.0) + c_p = JuMP.@constraint(wm.model, inv(L)*dhp <= sum(slopes_p .* qp_des_pipe)) + + dhn = var(wm, n, :dhn_des_pipe, a) + qn_des_pipe = var(wm, n, :qn_des_pipe, a) + qn_des_pipe_ub = JuMP.upper_bound.(qn_des_pipe) + slopes_n = pipe_resistances .* qn_des_pipe_ub.^(alpha - 1.0) + c_n = JuMP.@constraint(wm.model, inv(L)*dhn <= sum(slopes_n .* qn_des_pipe)) append!(con(wm, n, :head_loss)[a], [c_p, c_n]) end @@ -433,56 +435,82 @@ function constraint_resistance_selection_des(wm::AbstractDirectedModel, n::Int, for (r_id, r) in enumerate(pipe_resistances) x_res = var(wm, n, :x_res, a)[r_id] - qp_des = var(wm, n, :qp_des, a)[r_id] - qp_ub = JuMP.upper_bound(qp_des) - c_p = JuMP.@constraint(wm.model, qp_des <= qp_ub * x_res) + qp_des_pipe = var(wm, n, :qp_des_pipe, a)[r_id] + qp_ub = JuMP.upper_bound(qp_des_pipe) + c_p = JuMP.@constraint(wm.model, qp_des_pipe <= qp_ub * x_res) - qn_des = var(wm, n, :qn_des, a)[r_id] - qn_ub = JuMP.upper_bound(qn_des) - c_n = JuMP.@constraint(wm.model, qn_des <= qn_ub * x_res) + qn_des_pipe = var(wm, n, :qn_des_pipe, a)[r_id] + qn_ub = JuMP.upper_bound(qn_des_pipe) + c_n = JuMP.@constraint(wm.model, qn_des_pipe <= qn_ub * x_res) append!(con(wm, n, :head_loss)[a], [c_p, c_n]) end 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}, +function _gather_directionality_data( + wm::AbstractDirectedModel, n::Int, check_valve_fr::Array{Int64,1}, check_valve_to::Array{Int64,1}, pipe_fr::Array{Int64,1}, pipe_to::Array{Int64,1}, - 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}) + des_pipe_fr::Array{Int64,1}, des_pipe_to::Array{Int64,1}, pump_fr::Array{Int64,1}, + pump_to::Array{Int64,1}, pressure_reducing_valve_fr::Array{Int64,1}, + pressure_reducing_valve_to::Array{Int64,1}, shutoff_valve_fr::Array{Int64,1}, + shutoff_valve_to::Array{Int64,1}) # 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_pipe = var(wm, n, :y_pipe) + y_des_pipe = var(wm, n, :y_des_pipe) + y_pump = var(wm, n, :y_pump) y_pressure_reducing_valve = var(wm, n, :y_pressure_reducing_valve) y_shutoff_valve = var(wm, n, :y_shutoff_valve) - y_out = JuMP.@expression(wm.model, + sum_out = JuMP.@expression(wm.model, sum(y_check_valve[a] for a in check_valve_fr) + - sum(y_pipe[a] for a in pipe_fr) + sum(y_pump[a] for a in pump_fr) + + sum(y_pipe[a] for a in pipe_fr) + + sum(y_des_pipe[a] for a in des_pipe_fr) + + sum(y_pump[a] for a in pump_fr) + sum(y_pressure_reducing_valve[a] for a in pressure_reducing_valve_fr) + sum(y_shutoff_valve[a] for a in shutoff_valve_fr)) - y_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_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_pipe[a] for a in pipe_to) + + sum(y_des_pipe[a] for a in des_pipe_to) + + sum(y_pump[a] for a in pump_to) + sum(y_pressure_reducing_valve[a] for a in pressure_reducing_valve_to) + sum(y_shutoff_valve[a] for a in shutoff_valve_to)) - y_in_length = length(check_valve_to) + length(pipe_to) + length(pump_to) + - length(pressure_reducing_valve_to) + length(shutoff_valve_to) + # Get the in degree of node `i`. + in_length = length(check_valve_to) + length(pipe_to) + length(des_pipe_to) + + length(pump_to) + length(pressure_reducing_valve_to) + length(shutoff_valve_to) + + # Get the out degree of node `i`. + out_length = length(check_valve_fr) + length(pipe_fr) + length(des_pipe_fr) + + length(pump_fr) + length(pressure_reducing_valve_fr) + length(shutoff_valve_fr) + + return sum_in, sum_out, in_length, out_length +end + + +"Constraint to ensure at least one direction is set to take flow away from a source." +function constraint_node_directionality( + wm::AbstractDirectedModel, n::Int, i::Int, check_valve_fr::Array{Int64,1}, + check_valve_to::Array{Int64,1}, pipe_fr::Array{Int64,1}, pipe_to::Array{Int64,1}, + des_pipe_fr::Array{Int64,1}, des_pipe_to::Array{Int64,1}, pump_fr::Array{Int64,1}, + pump_to::Array{Int64,1}, pressure_reducing_valve_fr::Array{Int64,1}, + pressure_reducing_valve_to::Array{Int64,1}, shutoff_valve_fr::Array{Int64,1}, + shutoff_valve_to::Array{Int64,1}) + # Gather data required to build the constraint. + sum_in, sum_out, in_length, out_length = _gather_directionality_data( + wm, n, check_valve_fr, check_valve_to, pipe_fr, pipe_to, des_pipe_fr, des_pipe_to, + pump_fr, pump_to, pressure_reducing_valve_fr, pressure_reducing_valve_to, + shutoff_valve_fr, shutoff_valve_to) # 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) + if out_length == 1 && in_length == 1 + c = JuMP.@constraint(wm.model, sum_out - sum_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) + elseif in_length + out_length == 2 && in_length*out_length == 0 + c = JuMP.@constraint(wm.model, sum_out + sum_in == 1.0) con(wm, n, :node_directionality)[i] = c end end @@ -492,32 +520,18 @@ end function constraint_source_directionality( wm::AbstractDirectedModel, n::Int, i::Int, check_valve_fr::Array{Int64,1}, check_valve_to::Array{Int64,1}, pipe_fr::Array{Int64,1}, pipe_to::Array{Int64,1}, - 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_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) + des_pipe_fr::Array{Int64,1}, des_pipe_to::Array{Int64,1}, pump_fr::Array{Int64,1}, + pump_to::Array{Int64,1}, pressure_reducing_valve_fr::Array{Int64,1}, + pressure_reducing_valve_to::Array{Int64,1}, shutoff_valve_fr::Array{Int64,1}, + shutoff_valve_to::Array{Int64,1}) + # Gather data required to build the constraint. + sum_in, sum_out, in_length, out_length = _gather_directionality_data( + wm, n, check_valve_fr, check_valve_to, pipe_fr, pipe_to, des_pipe_fr, des_pipe_to, + pump_fr, pump_to, pressure_reducing_valve_fr, pressure_reducing_valve_to, + shutoff_valve_fr, shutoff_valve_to) # Add the source flow direction constraint. - c = JuMP.@constraint(wm.model, y_out - y_in >= 1.0 - y_in_length) + c = JuMP.@constraint(wm.model, sum_out - sum_in >= 1.0 - in_length) con(wm, n, :source_directionality)[i] = c end @@ -526,31 +540,17 @@ end function constraint_sink_directionality( wm::AbstractDirectedModel, n::Int, i::Int, check_valve_fr::Array{Int64,1}, check_valve_to::Array{Int64,1}, pipe_fr::Array{Int64,1}, pipe_to::Array{Int64,1}, - 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_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_out_length = length(check_valve_fr) + length(pipe_fr) + length(pump_fr) + - length(pressure_reducing_valve_fr) + length(shutoff_valve_fr) + des_pipe_fr::Array{Int64,1}, des_pipe_to::Array{Int64,1}, pump_fr::Array{Int64,1}, + pump_to::Array{Int64,1}, pressure_reducing_valve_fr::Array{Int64,1}, + pressure_reducing_valve_to::Array{Int64,1}, shutoff_valve_fr::Array{Int64,1}, + shutoff_valve_to::Array{Int64,1}) + # Gather data required to build the constraint. + sum_in, sum_out, in_length, out_length = _gather_directionality_data( + wm, n, check_valve_fr, check_valve_to, pipe_fr, pipe_to, des_pipe_fr, des_pipe_to, + pump_fr, pump_to, pressure_reducing_valve_fr, pressure_reducing_valve_to, + shutoff_valve_fr, shutoff_valve_to) # Add the sink flow direction constraint. - c = JuMP.@constraint(wm.model, y_in - y_out >= 1.0 - y_out_length) + c = JuMP.@constraint(wm.model, sum_in - sum_out >= 1.0 - out_length) con(wm, n, :sink_directionality)[i] = c end diff --git a/src/form/micp.jl b/src/form/micp.jl index c8545681..cc6389d7 100644 --- a/src/form/micp.jl +++ b/src/form/micp.jl @@ -4,11 +4,11 @@ function constraint_pipe_head_loss_des(wm::AbstractMICPModel, n::Int, a::Int, alpha::Float64, node_fr::Int, node_to::Int, L::Float64, pipe_resistances) # Collect head difference variables. - dhp, dhn = var(wm, n, :dhp_pipe, a), var(wm, n, :dhn_pipe, a) + dhp, dhn = var(wm, n, :dhp_des_pipe, a), var(wm, n, :dhn_des_pipe, a) for (r_id, r) in enumerate(pipe_resistances) # Collect directed flow variables. - qp, qn = var(wm, n, :qp_des, a)[r_id], var(wm, n, :qn_des, a)[r_id] + qp, qn = var(wm, n, :qp_des_pipe, a)[r_id], var(wm, n, :qn_des_pipe, a)[r_id] # Build the relaxed head loss constraints. c_p = JuMP.@NLconstraint(wm.model, r * head_loss(qp) <= inv(L) * dhp) diff --git a/src/form/milp.jl b/src/form/milp.jl index 69710642..68532f67 100644 --- a/src/form/milp.jl +++ b/src/form/milp.jl @@ -7,9 +7,9 @@ function variable_flow_piecewise_weights(wm::AbstractMILPModel; nw::Int=wm.cnw, if pipe_breakpoints > 0 # Create weights involved in convex combination constraints for pipes. lambda_pipe = var(wm, nw)[:lambda_pipe] = JuMP.@variable(wm.model, - [a in ids(wm, nw, :pipe_fixed), k in 1:pipe_breakpoints], + [a in ids(wm, nw, :pipe), k in 1:pipe_breakpoints], base_name="$(nw)_lambda", lower_bound=0.0, upper_bound=1.0, - start=comp_start_value(ref(wm, nw, :pipe_fixed, a), "lambda_start", k)) + start=comp_start_value(ref(wm, nw, :pipe, a), "lambda_start", k)) # Create weights involved in convex combination constraints for check valves. lambda_check_valve = var(wm, nw)[:lambda_check_valve] = JuMP.@variable(wm.model, @@ -22,6 +22,14 @@ function variable_flow_piecewise_weights(wm::AbstractMILPModel; nw::Int=wm.cnw, [a in ids(wm, nw, :shutoff_valve), k in 1:pipe_breakpoints], base_name="$(nw)_lambda", lower_bound=0.0, upper_bound=1.0, start=comp_start_value(ref(wm, nw, :shutoff_valve, a), "lambda_start", k)) + + # Create weights involved in convex combination constraints. + n_r = Dict(a=>length(ref(wm, nw, :resistance, a)) for a in ids(wm, nw, :des_pipe)) + + lambda_des_pipe = var(wm, nw)[:lambda_des_pipe] = JuMP.@variable(wm.model, + [a in ids(wm, nw, :des_pipe), r in 1:n_r[a], k in 1:pipe_breakpoints], + base_name="$(nw)_lambda", lower_bound=0.0, upper_bound=1.0, + start=comp_start_value(ref(wm, nw, :des_pipe, a), "lambda_start", r)) end pump_breakpoints = get(wm.ext, :pump_breakpoints, 0) @@ -36,21 +44,6 @@ function variable_flow_piecewise_weights(wm::AbstractMILPModel; nw::Int=wm.cnw, end -function variable_flow_piecewise_weights_des(wm::AbstractMILPModel; nw::Int=wm.cnw, report::Bool=false) - pipe_breakpoints = get(wm.ext, :pipe_breakpoints, 0) - - if pipe_breakpoints > 0 - n_r = Dict(a=>length(ref(wm, nw, :resistance, a)) for a in ids(wm, nw, :pipe_des)) - - # Create weights involved in convex combination constraints. - lambda_pipe = var(wm, nw)[:lambda_pipe] = JuMP.@variable(wm.model, - [a in ids(wm, nw, :pipe_des), r in 1:n_r[a], k in 1:pipe_breakpoints], - base_name="$(nw)_lambda", lower_bound=0.0, upper_bound=1.0, - start=comp_start_value(ref(wm, nw, :pipe_des, a), "lambda_start", r)) - end -end - - function variable_flow_piecewise_adjacency(wm::AbstractMILPModel; nw::Int=wm.cnw, report::Bool=false) pipe_breakpoints = get(wm.ext, :pipe_breakpoints, 0) @@ -69,6 +62,11 @@ function variable_flow_piecewise_adjacency(wm::AbstractMILPModel; nw::Int=wm.cnw x_pw_shutoff_valve = var(wm, nw)[:x_pw_shutoff_valve] = JuMP.@variable(wm.model, [a in ids(wm, nw, :shutoff_valve), k in 1:pipe_breakpoints-1], base_name="$(nw)_x_pw", binary=true, start=comp_start_value(ref(wm, nw, :shutoff_valve, a), "x_pw_start")) + + # Create binary variables for design pipe convex combination constraints. + x_pw_des_pipe = var(wm, nw)[:x_pw_des_pipe] = JuMP.@variable(wm.model, + [a in ids(wm, nw, :des_pipe), k in 1:pipe_breakpoints-1], base_name="$(nw)_x_pw", binary=true, + start=comp_start_value(ref(wm, nw, :des_pipe, a), "x_pw_start")) end pump_breakpoints = get(wm.ext, :pump_breakpoints, 0) @@ -92,14 +90,16 @@ function variable_flow(wm::AbstractMILPModel; nw::Int=wm.cnw, bounded::Bool=true # Create variables required for convex combination piecewise approximation. variable_flow_piecewise_weights(wm, nw=nw) variable_flow_piecewise_adjacency(wm, nw=nw) + + # Create common flow variables. + variable_flow_des_common(wm, nw=nw, bounded=bounded, report=report) end -"Creates network design flow variables for `MILP` formulations (`q_des`, `lambda`, `x_pw`)." +"Creates network design flow variables for `MILP` formulations (`q_des_pipe`, `lambda`, `x_pw`)." function variable_flow_des(wm::AbstractMILPModel; nw::Int=wm.cnw, bounded::Bool=true, report::Bool=true) # Create common flow variables. variable_flow_des_common(wm, nw=nw, bounded=bounded, report=report) - variable_flow_piecewise_weights_des(wm, nw=nw) end @@ -241,14 +241,14 @@ function constraint_pipe_head_loss_des(wm::AbstractMILPModel, n::Int, a::Int, al # Get common variables and data. n_b = pipe_breakpoints - lambda, x_pw = var(wm, n, :lambda_pipe), var(wm, n, :x_pw_pipe) + lambda, x_pw = var(wm, n, :lambda_des_pipe), var(wm, n, :x_pw_des_pipe) h_i, h_j = var(wm, n, :h, node_fr), var(wm, n, :h, node_to) # Initialize flow expression in the head loss relationship. q_loss_expr = JuMP.AffExpr(0.0) for (r_id, r) in enumerate(resistances) - q, x_res = var(wm, n, :q_des, a)[r_id], var(wm, n, :x_res, a)[r_id] + q, x_res = var(wm, n, :q_des_pipe, a)[r_id], var(wm, n, :x_res, a)[r_id] q_lb, q_ub = JuMP.lower_bound(q), JuMP.upper_bound(q) # Add the first required SOS constraints. diff --git a/src/form/milpr.jl b/src/form/milpr.jl index 254737b8..981e067d 100644 --- a/src/form/milpr.jl +++ b/src/form/milpr.jl @@ -72,7 +72,7 @@ function constraint_pipe_head_loss_des(wm::AbstractMILPRModel, n::Int, a::Int, a for (r_id, r) in enumerate(resistances) # Add constraints corresponding to positive outer approximations. - qp, dhp = var(wm, n, :qp_des, a)[r_id], var(wm, n, :dhp_pipe, a) + qp, dhp = var(wm, n, :qp_des_pipe, a)[r_id], var(wm, n, :dhp_des_pipe, a) qp_ub = JuMP.upper_bound(qp) for qp_hat in range(0.0, stop=qp_ub, length=pipe_breakpoints) @@ -82,7 +82,7 @@ function constraint_pipe_head_loss_des(wm::AbstractMILPRModel, n::Int, a::Int, a end # Add constraints corresponding to negative outer approximations. - qn, dhn = var(wm, n, :qn_des, a)[r_id], var(wm, n, :dhn_pipe, a) + qn, dhn = var(wm, n, :qn_des_pipe, a)[r_id], var(wm, n, :dhn_des_pipe, a) qn_ub = JuMP.upper_bound(qn) for qn_hat in range(0.0, stop=qn_ub, length=pipe_breakpoints) diff --git a/src/form/nlp.jl b/src/form/nlp.jl index 27e0088f..b1f271f6 100644 --- a/src/form/nlp.jl +++ b/src/form/nlp.jl @@ -12,7 +12,7 @@ #end -"Creates network design flow variables for `NLP` formulations (`q_des`, `x_des`)." +"Creates network design flow variables for `NLP` formulations (`q_des_pipe`, `x_des`)." function variable_flow_des(wm::AbstractNLPModel; nw::Int=wm.cnw, bounded::Bool=true, report::Bool=true) variable_flow_des_common(wm, nw=nw, bounded=bounded, report=report) end @@ -78,7 +78,7 @@ end "Add head loss constraints for design pipes (without check valves) in `NLP` formulations." function constraint_pipe_head_loss_des(wm::AbstractNLPModel, n::Int, a::Int, alpha::Float64, node_fr::Int, node_to::Int, L::Float64, res) # Gather common flow and head variables, as well as design indices. - q, R = var(wm, n, :q_des, a), 1:length(res) + q, R = var(wm, n, :q_des_pipe, a), 1:length(res) h_i, h_j = var(wm, n, :h, node_fr), var(wm, n, :h, node_to) # Add the nonconvex, design-expanded head loss constraint. diff --git a/src/form/undirected_flow.jl b/src/form/undirected_flow.jl index 68d66cfc..de383303 100644 --- a/src/form/undirected_flow.jl +++ b/src/form/undirected_flow.jl @@ -9,6 +9,9 @@ function variable_flow(wm::AbstractUndirectedModel; nw::Int=wm.cnw, bounded::Boo # Create directed flow (`qp` and `qn`) variables for each component. variable_component_flow(wm, name; nw=nw, bounded=bounded, report=report) end + + # Create flow-related variables for design components. + variable_flow_des_common(wm, nw=nw, bounded=bounded, report=report) end @@ -40,33 +43,33 @@ end "Create common network design flow variables for undirected flow formulations." function variable_flow_des_common(wm::AbstractUndirectedModel; nw::Int=wm.cnw, bounded::Bool=true, report::Bool=true) - # Create dictionary for undirected design flow variables (i.e., q_des). - q_des = var(wm, nw)[:q_des] = Dict{Int,Array{JuMP.VariableRef}}() + # Create dictionary for undirected design flow variables (i.e., q_des_pipe). + q_des_pipe = var(wm, nw)[:q_des_pipe] = Dict{Int,Array{JuMP.VariableRef}}() # Initialize the variables. (The default start value of _q_eps is crucial.) - for a in ids(wm, nw, :pipe_des) - var(wm, nw, :q_des)[a] = JuMP.@variable(wm.model, - [r in 1:length(ref(wm, nw, :resistance, a))], base_name="$(nw)_q_des", - start=comp_start_value(ref(wm, nw, :pipe_des, a), "q_des_start", r, _q_eps)) + for a in ids(wm, nw, :des_pipe) + var(wm, nw, :q_des_pipe)[a] = JuMP.@variable(wm.model, + [r in 1:length(ref(wm, nw, :resistance, a))], base_name="$(nw)_q_des_pipe", + start=comp_start_value(ref(wm, nw, :des_pipe, a), "q_des_pipe_start", r, _q_eps)) end if bounded # If the variables are bounded, apply the bounds. q_lb, q_ub = calc_flow_bounds(wm, nw) - for a in ids(wm, nw, :pipe_des) + for a in ids(wm, nw, :des_pipe) for r in 1:length(ref(wm, nw, :resistance, a)) - JuMP.set_lower_bound(q_des[a][r], q_lb["pipe"][a][r]) - JuMP.set_upper_bound(q_des[a][r], q_ub["pipe"][a][r]) + JuMP.set_lower_bound(q_des_pipe[a][r], q_lb["des_pipe"][a][r]) + JuMP.set_upper_bound(q_des_pipe[a][r], q_ub["des_pipe"][a][r]) end end end - # Create expressions capturing the relationships among q and q_des. - q = var(wm, nw)[:q] = JuMP.@expression( - wm.model, [a in ids(wm, nw, :pipe_des)], sum(var(wm, nw, :q_des, a))) + # Create expressions capturing the relationships among q and q_des_pipe. + q = var(wm, nw)[:q_des_pipe_sum] = JuMP.@expression( + wm.model, [a in ids(wm, nw, :des_pipe)], sum(var(wm, nw, :q_des_pipe, a))) # Initialize the solution reporting data structures. - report && sol_component_value(wm, nw, :pipe_des, :q, ids(wm, nw, :pipe_des), q) + report && sol_component_value(wm, nw, :des_pipe, :q, ids(wm, nw, :des_pipe), q) # Create resistance binary variables. variable_resistance(wm, nw=nw) @@ -79,14 +82,14 @@ function constraint_resistance_selection_des(wm::AbstractUndirectedModel, n::Int append!(con(wm, n, :head_loss)[a], [c]) for r in 1:length(pipe_resistances) - q_des = var(wm, n, :q_des, a)[r] + q_des_pipe = var(wm, n, :q_des_pipe, a)[r] x_res = var(wm, n, :x_res, a)[r] - q_des_lb = JuMP.lower_bound(q_des) - c_lb = JuMP.@constraint(wm.model, q_des >= q_des_lb * x_res) + q_des_pipe_lb = JuMP.lower_bound(q_des_pipe) + c_lb = JuMP.@constraint(wm.model, q_des_pipe >= q_des_pipe_lb * x_res) - q_des_ub = JuMP.upper_bound(q_des) - c_ub = JuMP.@constraint(wm.model, q_des <= q_des_ub * x_res) + q_des_pipe_ub = JuMP.upper_bound(q_des_pipe) + c_ub = JuMP.@constraint(wm.model, q_des_pipe <= q_des_pipe_ub * x_res) append!(con(wm, n, :head_loss)[a], [c_lb, c_ub]) end @@ -180,30 +183,36 @@ 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}) + wm::AbstractUndirectedModel, n::Int, i::Int, check_valve_fr::Array{Int64,1}, + check_valve_to::Array{Int64,1}, pipe_fr::Array{Int64,1}, pipe_to::Array{Int64,1}, + des_pipe_fr::Array{Int64,1}, des_pipe_to::Array{Int64,1}, pump_fr::Array{Int64,1}, + pump_to::Array{Int64,1}, pressure_reducing_valve_fr::Array{Int64,1}, + pressure_reducing_valve_to::Array{Int64,1}, shutoff_valve_fr::Array{Int64,1}, + shutoff_valve_to::Array{Int64,1}) # 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}, - 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}) + wm::AbstractUndirectedModel, n::Int, i::Int, check_valve_fr::Array{Int64,1}, + check_valve_to::Array{Int64,1}, pipe_fr::Array{Int64,1}, pipe_to::Array{Int64,1}, + des_pipe_fr::Array{Int64,1}, des_pipe_to::Array{Int64,1}, pump_fr::Array{Int64,1}, + pump_to::Array{Int64,1}, pressure_reducing_valve_fr::Array{Int64,1}, + pressure_reducing_valve_to::Array{Int64,1}, shutoff_valve_fr::Array{Int64,1}, + shutoff_valve_to::Array{Int64,1}) # For undirected formulations, there are no constraints, here. end + function constraint_source_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}) + wm::AbstractUndirectedModel, n::Int, i::Int, check_valve_fr::Array{Int64,1}, + check_valve_to::Array{Int64,1}, pipe_fr::Array{Int64,1}, pipe_to::Array{Int64,1}, + des_pipe_fr::Array{Int64,1}, des_pipe_to::Array{Int64,1}, pump_fr::Array{Int64,1}, + pump_to::Array{Int64,1}, pressure_reducing_valve_fr::Array{Int64,1}, + pressure_reducing_valve_to::Array{Int64,1}, shutoff_valve_fr::Array{Int64,1}, + shutoff_valve_to::Array{Int64,1}) # For undirected formulations, there are no constraints, here. end diff --git a/src/prob/des.jl b/src/prob/des.jl index 676296d9..41504545 100644 --- a/src/prob/des.jl +++ b/src/prob/des.jl @@ -10,7 +10,6 @@ function build_des(wm::AbstractWaterModel) variable_head(wm) variable_head_gain(wm) variable_flow(wm) - variable_flow_des(wm) # Component-specific variables. variable_reservoir(wm) @@ -18,29 +17,32 @@ function build_des(wm::AbstractWaterModel) # Add the network design objective. objective_des(wm) - for (a, pipe) in ref(wm, :pipe_fixed) + # 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 fixed (non-design) pipes. + for (a, pipe) in ref(wm, :pipe) constraint_pipe_head_loss(wm, a) end + # Head loss along design pipes. + for (a, pipe) in ref(wm, :des_pipe) + constraint_pipe_head_loss_des(wm, a) + end + + # Head loss along pipes with check valves. for (a, check_valve) in ref(wm, :check_valve) constraint_check_valve_head_loss(wm, a) end + # Head loss along pipes with shutoff valves. for (a, shutoff_valve) in ref(wm, :shutoff_valve) constraint_shutoff_valve_head_loss(wm, a) end - # Head loss along design pipes. - for (a, pipe) in ref(wm, :pipe_des) - constraint_pipe_head_loss_des(wm, a) - end - - # 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. for (i, reservoir) in ref(wm, :reservoir) constraint_reservoir_head(wm, reservoir["node"]) @@ -49,9 +51,13 @@ function build_des(wm::AbstractWaterModel) # Constrain flow directions based on demand. for (i, junction) in ref(wm, :junction) - if junction["demand"] > 0.0 + if !junction["dispatchable"] && junction["demand"] > 0.0 constraint_sink_directionality(wm, junction["node"]) - elseif junction["demand"] < 0.0 + elseif junction["dispatchable"] && junction["demand_min"] > 0.0 + constraint_sink_directionality(wm, junction["node"]) + elseif !junction["dispatchable"] && junction["demand"] < 0.0 + constraint_source_directionality(wm, junction["node"]) + elseif junction["dispatchable"] && junction["demand_max"] < 0.0 constraint_source_directionality(wm, junction["node"]) end end diff --git a/test/des.jl b/test/des.jl index 8e91d805..a47031e6 100644 --- a/test/des.jl +++ b/test/des.jl @@ -1,5 +1,5 @@ @testset "Network Design Problems (Single Network)" begin - network = parse_file("../test/data/epanet/shamir.inp") + network = parse_file("../examples/data/epanet/shamir.inp") modifications = parse_file("../test/data/json/shamir-reduced.json") _IM.update_data!(network, modifications) diff --git a/test/io.jl b/test/io.jl index 36a0c14e..0dc49c95 100644 --- a/test/io.jl +++ b/test/io.jl @@ -14,8 +14,8 @@ klmod_name = "Global Water Full network - Peak Day (Avg * 1.9)" @test klmod_data["name"] == klmod_name - shamir_data = WaterModels.parse_file("../test/data/epanet/shamir.inp") - shamir_name = "shamir -- Bragalli, D'Ambrosio, Lee, Lodi, Toth (2008)" + shamir_data = WaterModels.parse_file("../examples/data/epanet/shamir.inp") + shamir_name = "Shamir (Two-loop) Water Network" @test shamir_data["name"] == shamir_name shamir_ts_data = WaterModels.parse_file("../test/data/epanet/shamir-ts.inp") @@ -24,18 +24,18 @@ end @testset "parse_file (.json)" begin - data = WaterModels.parse_file("../test/data/json/shamir.json") + data = WaterModels.parse_file("../examples/data/json/shamir.json") pipe_1_max_velocity = data["pipe"]["1"]["maximumVelocity"] @test pipe_1_max_velocity == 2.0 end @testset "parse_file (invalid extension)" begin - path = "../test/data/json/shamir.data" + path = "../examples/data/json/shamir.data" @test_throws ErrorException WaterModels.parse_file(path) end @testset "parse_json" begin - data = WaterModels.parse_json("../test/data/json/shamir.json") + data = WaterModels.parse_json("../examples/data/json/shamir.json") pipe_1_max_velocity = data["pipe"]["1"]["maximumVelocity"] @test pipe_1_max_velocity == 2.0 end diff --git a/test/warm_start.jl b/test/warm_start.jl index 64e4418f..e925062b 100644 --- a/test/warm_start.jl +++ b/test/warm_start.jl @@ -1,6 +1,6 @@ @testset "Warm Start Functionality" begin @testset "Shamir network, wf, NLP formulation." begin - data = WaterModels.parse_file("../test/data/epanet/shamir.inp") + data = WaterModels.parse_file("../examples/data/epanet/shamir.inp") cold_result = WaterModels.run_wf(data, NLPWaterModel, ipopt) @test cold_result["termination_status"] == LOCALLY_SOLVED _IM.update_data!(data, cold_result["solution"]) From 47bcb56773eb4472501959cbd48cb22301b2c203 Mon Sep 17 00:00:00 2001 From: Byron Tasseff Date: Sat, 22 Aug 2020 12:57:35 -0600 Subject: [PATCH 2/3] Minor modification changes. --- docs/src/math-model.md | 115 +++++++------------------------------ docs/src/quickguide.md | 8 +-- docs/src/specifications.md | 82 +++++++++++++------------- 3 files changed, 64 insertions(+), 141 deletions(-) diff --git a/docs/src/math-model.md b/docs/src/math-model.md index 882715bd..296c3092 100644 --- a/docs/src/math-model.md +++ b/docs/src/math-model.md @@ -1,26 +1,20 @@ # Mathematical Models in WaterModels ## Notation for Sets -A water distribution network is represented by a directed graph $\mathcal{G} := (\mathcal{N}, \mathcal{L})$, where $\mathcal{N}$ is the set of nodes (e.g., [junctions](http://wateranalytics.org/EPANET/_juncs_page.html) and [reservoirs](http://wateranalytics.org/EPANET/_resv_page.html)) and $\mathcal{L}$ is the set of arcs (conventionally "links," e.g., [pipes](http://wateranalytics.org/EPANET/_pipes_page.html) and [valves](http://wateranalytics.org/EPANET/_valves_page.html)). +A water distribution network is represented by a directed graph $\mathcal{G} := (\mathcal{N}, \mathcal{L})$, where $\mathcal{N}$ is the set of nodes and $\mathcal{L}$ is the set of arcs (conventionally "links," e.g., [pipes](http://wateranalytics.org/EPANET/_pipes_page.html) and [valves](http://wateranalytics.org/EPANET/_valves_page.html)). Temporal evolution of the network is represented by a set $\mathcal{K}$, denoting the set of all time steps considered. -The set of junctions is denoted by $\mathcal{J} \subset \mathcal{N}$, the set of reservoirs (or sources) by $\mathcal{R} \subset \mathcal{N}$, and the set of tanks by $\mathcal{T} \subset \mathcal{N}$. -The set of pipes in the network is denoted by $\mathcal{A} \subset \mathcal{L}$, the set of pipes with check valves by $\mathcal{C} \subset \mathcal{A}$, and the set of pumps by $\mathcal{P} \subset \mathcal{L}$. -The set of links incident to node $i \in \mathcal{N}$, where $i$ is the tail of the arc, is denoted by $\delta^{+}_{i} := \{(i, j) \in \mathcal{L}\}$. -The set of links incident to node $i \in \mathcal{N}$, where $i$ is the head of the arc, is denoted by $\delta^{-}_{i} := \{(j, i) \in \mathcal{L}\}$. In summary, the following sets are commonly used when defining a WaterModels problem formulation: | Notation | WaterModels Translation | Description | | :-------------------------------------- | :----------------------------- | :------------------------- | -| $\mathcal{N}$ | `wm.ref[:nw][n][:node]` | nodes | +| $\mathcal{N}$ | `wm.ref[:nw][n][:node]` | nodes (to which nodal-type components are attached) | | $\mathcal{K}$ | `nw_ids(wm)` | time indices (multinetwork indices labeled by `n`) | -| $\mathcal{J} \subset \mathcal{N}$ | `wm.ref[:nw][n][:junction]` | [junctions](http://wateranalytics.org/EPANET/_juncs_page.html) | -| $\mathcal{R} \subset \mathcal{N}$ | `wm.ref[:nw][n][:reservoir]` | [reservoirs](http://wateranalytics.org/EPANET/_resv_page.html) | -| $\mathcal{T} \subset \mathcal{N}$ | `wm.ref[:nw][n][:tank]` | [tanks](http://wateranalytics.org/EPANET/_tanks_page.html) | +| $\mathcal{J}$ | `wm.ref[:nw][n][:junction]` | [junctions](http://wateranalytics.org/EPANET/_juncs_page.html) | +| $\mathcal{R}$ | `wm.ref[:nw][n][:reservoir]` | [reservoirs](http://wateranalytics.org/EPANET/_resv_page.html) | +| $\mathcal{T}$ | `wm.ref[:nw][n][:tank]` | [tanks](http://wateranalytics.org/EPANET/_tanks_page.html) | | $\mathcal{A} \subset \mathcal{L}$ | `wm.ref[:nw][n][:pipe]` | [pipes](http://wateranalytics.org/EPANET/_pipes_page.html) | | $\mathcal{C} \subset \mathcal{L}$ | `wm.ref[:nw][n][:check_valve]` | [pipes with check valves](http://wateranalytics.org/EPANET/_pipes_page.html) | | $\mathcal{P} \subset \mathcal{L}$ | `wm.ref[:nw][n][:pump]` | [pumps](http://wateranalytics.org/EPANET/_pumps_page.html) | -| $\delta_{i}^{+} \subset \mathcal{A}$ | `wm.ref[:nw][n][:link_to][i]` | links "from" node $i$ | -| $\delta_{i}^{-} \subset \mathcal{A}$ | `wm.ref[:nw][n][:link_fr][i]` | links "to" node $i$ | ## Physical Feasibility ### Nodes @@ -37,10 +31,14 @@ In summary, the following sets are commonly used when defining a WaterModels pro #### Pipes with Check Valves +#### Pipes with Shutoff Valves + +#### Pressure Reducing Valves + #### Pumps ### Satisfaction of Flow Bounds -For each arc $(i, j) \in \mathcal{A}$, a variable $q_{ij}$ is used to represent the volumetric flow of water across the arc (in $\textrm{m}^{3}/\textrm{s}$). +For each arc $(i, j) \in \mathcal{L}$, a variable $q_{ij}$ is used to represent the volumetric flow of water across the arc (in $\textrm{m}^{3}/\textrm{s}$). When $q_{ij}$ is positive, flow on arc $(i, j)$ travels from node $i$ to node $j$. When $q_{ij}$ is negative, flow travels from node $j$ to node $i$. The absolute value of flow along the arc can be bounded by physical capacity, engineering judgment, or network analysis. @@ -54,44 +52,28 @@ where $D_{ij}$ is the diameter of pipe $(i, j)$ and $v^{\max}_{ij}$ is the maxim ### Satisfaction of Head Bounds Each node potential is denoted by $h_{i}$, $i \in \mathcal{N}$, and represents the hydraulic head in units of length ($\textrm{m}$). The hydraulic head assimilates the elevation and pressure heads at each node, while the velocity head can typically be neglected. -For each reservoir $i \in \mathcal{S}$, the hydraulic head is assumed to be fixed at a value $h_{i}^{\textrm{src}}$, i.e., +For each reservoir $i \in \mathcal{R}$, the hydraulic head is assumed to be fixed at a value $h_{i}^{\textrm{src}}$, i.e., ```math - h_{i} = h_{i}^{\textrm{src}}, \; \forall i \in \mathcal{S}. + h_{i} = h_{i}^{\textrm{src}}, \; \forall i \in \mathcal{R}. ``` For each junction $i \in \mathcal{J}$, a minimum hydraulic head $\underline{h}_{i}$, determined a priori, must first be satisfied. In the interest of tightening the optimization formulation, upper bounds on hydraulic heads can also typically be implied from other network data, e.g., ```math - \underline{h}_{i} \leq h_{i} \leq \overline{h}_{i} = \max_{i \in \mathcal{S}}\{h_{i}^{\textrm{src}}\}. + \underline{h}_{i} \leq h_{i} \leq \overline{h}_{i} = \max_{i \in \mathcal{R}}\{h_{i}^{\textrm{src}}\}. ``` -### Conservation of Flow at Non-supply Nodes -Flow must be delivered throughout the network to satisfy fixed demand, $q_{i}^{\textrm{dem}}$, at non-supply nodes, i.e., -```math - \sum_{(j, i) \in \mathcal{A}^{-}(i)} q_{ji} - \sum_{(i, j) \in \mathcal{A}^{+}(i)} q_{ij} = q_{i}^{\textrm{dem}}, ~ \forall i \in \mathcal{J}, -``` -where $\mathcal{A}^{-}(i)$ and $\mathcal{A}^{+}(i)$ are the sets of incoming and outgoing arcs of node $i$, respectively. - -### Conservation of Flow at Supply Nodes -The _outflow_ from each reservoir will be nonnegative by definition, i.e., -```math - \sum_{(i, j) \in \mathcal{A}^{+}(i)} q_{ij} - \sum_{(j, i) \in \mathcal{A}^{-}(i)} q_{ji} \geq 0, ~ \forall i \in \mathcal{S}. -``` -Additionally, an upper bound on the amount of flow delivered by a reservoir may be written -```math - \sum_{(i, j) \in \mathcal{A}^{+}(i)} q_{ij} - \sum_{(j, i) \in \mathcal{A}^{-}(i)} q_{ji} \leq \sum_{k \in \mathcal{J}} q^{\textrm{dem}}_{k}, ~ \forall i \in \mathcal{R}, -``` -i.e., a reservoir will never send more flow than the amount required to serve all demand. +### Conservation of Flow ### Head Loss Relationships -In water distribution networks, flow along an arc is induced by the difference in potential (head) between the two nodes that connect that arc. +In water distribution networks, flow along a pipe is induced by the difference in potential (head) between the two nodes that connect that pipe. The relationships that link flow and hydraulic head are commonly referred to as the "head loss equations" or "potential-flow constraints," and are generally of the form ```math h_{i} - h_{j} = \Phi_{ij}(q_{ij}), ``` where $\Phi_{ij} : \mathbb{R} \to \mathbb{R}$ is a strictly increasing function with rotational symmetry about the origin. -Embedding the above equation in a mathematical program clearly introduces non-convexity. -(That is, the function $\Phi_{ij}(q_{ij})$ is non-convex _and_ the relationship must be satisfied with equality.) -As such, different formulations primarily aim to effectively deal with these types of non-convex constraints in an optimization setting. +Embedding the above equation in a mathematical program clearly introduces nonconvexity. +(That is, the function $\Phi_{ij}(q_{ij})$ is nonconvex _and_ the relationship must be satisfied with equality.) +As such, different formulations primarily aim to effectively deal with these types of nonconvex constraints in an optimization setting. Explicit forms of the head loss equation include the [Darcy-Weisbach](https://en.wikipedia.org/wiki/Darcy-Weisbach_equation) equation, i.e., ```math @@ -129,8 +111,8 @@ and the Hazen-Williams resistance per unit length is r_{ij} = \frac{10.67}{\kappa_{ij}^{1.852} D_{ij}^{4.87}}. ``` -## Non-convex Nonlinear Program -The full non-convex formulation of the physical feasibility problem (NCNLP), which incorporates all requirements from [Physical Feasibility](#Physical-Feasibility-1), may be written as a system that satisfies the following constraints: +## Nonconvex Nonlinear Program +The full nonconvex formulation of the physical feasibility problem (NCNLP), which incorporates all requirements from [Physical Feasibility](#Physical-Feasibility-1), may be written as a system that satisfies the following constraints: ```math \begin{align} h_{i} - h_{j} &= L_{ij} r_{ij} q_{ij} \lvert q_{ij} \rvert^{\alpha}, ~ \forall (i, j) \in \mathcal{A} \label{eqn:ncnlp-head-loss} \\ @@ -140,62 +122,9 @@ The full non-convex formulation of the physical feasibility problem (NCNLP), whi \underline{q}_{ij} \leq q_{ij} &\leq \overline{q}_{ij}, ~ \forall (i, j) \in \mathcal{A} \label{eqn:ncnlp-flow-bounds}. \end{align} ``` -Here, Constraints $\eqref{eqn:ncnlp-head-loss}$ are [head loss relationships](#Head-Loss-Relationships-1), Constraints $\eqref{eqn:ncnlp-head-source}$ are [head bounds](#Satisfaction-of-Head-Bounds-1) at source nodes, Constraints $\eqref{eqn:ncnlp-flow-conservation}$ are [flow conservation constraints](#Conservation-of-Flow-at-Non-supply-Nodes-1), Constraints $\eqref{eqn:ncnlp-head-bounds}$ [head bounds](#Satisfaction-of-Head-Bounds-1) at junctions, and Constraints $\eqref{eqn:ncnlp-flow-bounds}$ are [flow bounds](#Satisfaction-of-Flow-Bounds-1). -Note that the sources of non-convexity and nonlinearity are Constraints $\eqref{eqn:ncnlp-head-loss}$. - -## Convex Nonlinear Program -Note that the sources of non-convexity and nonlinearity in the full [non-convex formulation of the physical feasibility problem](#Non-convex-Nonlinear-Program-1) are Constraints $\eqref{eqn:ncnlp-head-loss}$. -Because of the symmetry of the head loss relationship, the problem can be modeled instead as a disjunctive program. -Here, the disjunction arises from the direction of flow, i.e., -```math -\begin{equation} - \left[ - \begin{aligned}[c] - h_{i} - h_{j} &= L_{ij} r_{ij} q_{ij}^{1 + \alpha} \\ - q_{ij} &\geq 0 - \end{aligned} - \right] - \lor - \left[ - \begin{aligned}[c] - h_{i} - h_{j} &= L_{ij} r_{ij} (-q_{ij})^{1 + \alpha} \\ - q_{ij} &< 0 - \end{aligned} - \right], ~ \forall (i, j) \in \mathcal{A}, \label{eqn:dnlp-head-loss} -\end{equation} -``` -which replaces Constraints $\eqref{eqn:ncnlp-head-loss}$ in the [NCNLP formulation](#Non-convex-Nonlinear-Program-1). -To model the disjunction, each flow variable $q_{ij}$ can be decomposed into two nonnegative flow variables, $q_{ij}^{+}$ and $q_{ij}^{-}$, where $q_{ij} := q_{ij}^{+} - q_{ij}^{-}$. -With this in mind, the following _convex_ nonlinear program (CNLP) can be formulated, which is adapted from Section 3 of [_Global Optimization of Nonlinear Network Design_ by Raghunathan (2013)](https://epubs.siam.org/doi/abs/10.1137/110827387): -```math -\begin{align} - & \text{minimize} - & & \sum_{(i, j) \in \mathcal{A}} \frac{L_{ij} r_{ij}}{2 + \alpha} \left[(q_{ij}^{+})^{2 + \alpha} + (q_{ij}^{-})^{2 + \alpha}\right] - \sum_{i \in \mathcal{S}} h_{i}^{\textrm{src}} \left(\sum_{(i, j) \in \mathcal{A}^{-}(i)} (q_{ij}^{+} - q_{ij}^{-}) - \sum_{(j, i) \in \mathcal{A}^{+}(i)} (q_{ji}^{+} - q_{ji}^{-})\right) \\ - & \text{subject to} - & & \sum_{(j, i) \in \mathcal{A}^{-}(i)} (q_{ji}^{+} - q_{ji}^{-}) - \sum_{(i, j) \in \mathcal{A}^{+}(i)} (q_{ij}^{+} - q_{ij}^{-}) = q_{i}^{\textrm{dem}}, ~ \forall i \in \mathcal{J} \label{eqn:cnlp-flow-conservation} \\ - & & & q_{ij}^{+}, q_{ij}^{-} \geq 0, ~ \forall (i, j) \in \mathcal{A} \label{eqn:cnlp-flow-bounds}. -\end{align} -``` -$\renewcommand{\hat}[1]{\widehat{#1}}$ -Suppose that $\hat{\mathbf{q}}^{+}, \hat{\mathbf{q}}^{-} \in \mathbb{R}^{\lvert A \rvert}$ solves (CNLP) with the associated dual solution $\hat{\mathbf{h}} \in \mathbb{R}^{\lvert \mathcal{J} \rvert}$, corresponding to the flow conservation Constraints $\eqref{eqn:cnlp-flow-conservation}$, and $\hat{\mathbf{u}}^{+}, \hat{\mathbf{u}}^{-} \in \mathbb{R}^{\lvert \mathcal{A} \rvert}$, corresponding to the nonnegativity Constraints $\eqref{eqn:cnlp-flow-bounds}$. -This solution must satisfy the first-order necessary conditions -```math -\begin{align} - \hat{h}_{i} - \hat{h}_{j} &= L_{ij} r_{ij} \hat{q}_{ij} \lvert q_{ij} \rvert^{\alpha}, ~ \forall (i, j) \in \mathcal{A} \\ - h_{i} &= h_{i}^{\textrm{src}}, ~ \forall i \in \mathcal{S} \\ - \sum_{(j, i) \in \mathcal{A}^{-}(i)} q_{ji} - \sum_{(i, j) \in \mathcal{A}^{+}(i)} q_{ij} &= q_{i}^{\textrm{dem}}, ~ \forall i \in \mathcal{J} \\ - \underline{h}_{i} \leq h_{i} &\leq \overline{h}_{i}, ~ \forall i \in \mathcal{J} \\ - \underline{q}_{ij} \leq q_{ij} &\leq \overline{q}_{ij}, ~ \forall (i, j) \in \mathcal{A}. -\end{align} -``` +Here, Constraints $\eqref{eqn:ncnlp-head-loss}$ are [head loss relationships](#Head-Loss-Relationships-1), Constraints $\eqref{eqn:ncnlp-head-source}$ are [head bounds](#Satisfaction-of-Head-Bounds-1) at source nodes, Constraints $\eqref{eqn:ncnlp-flow-conservation}$ are [flow conservation constraints](#Conservation-of-Flow), Constraints $\eqref{eqn:ncnlp-head-bounds}$ [head bounds](#Satisfaction-of-Head-Bounds-1) at junctions, and Constraints $\eqref{eqn:ncnlp-flow-bounds}$ are [flow bounds](#Satisfaction-of-Flow-Bounds-1). +Note that the sources of nonconvexity and nonlinearity are Constraints $\eqref{eqn:ncnlp-head-loss}$. ## Mixed-integer Convex Program ## Mixed-integer Linear Program - -## Network Design -Currently, the primary formulation focuses on the problem of optimally designing a water distribution network. -More specifically, given a network consisting of reservoirs, junctions, and pipes, the problem aims to select the most cost-effecient resistance from a discrete set of resistances for each pipe to meet demand over the entire network. -The set of all possible resistances for a given pipe $(i, j) \in \mathcal{A}$ is denoted by $\mathcal{R}_{ij}$, where each resistance is denoted by $r \in \mathcal{R}_{ij}$. -A binary variable $x^{r}_{ijr}$ is associated with each of these diameters to model the decision, i.e., $x_{ijr}^{r} = 1$ if $r$ is selected to serve as the pipe resistance, and $x_{ijr}^{r} = 0$ otherwise. -The cost per unit length of installing a pipe of resistance $r$ is denoted by $c_{ijr}$. diff --git a/docs/src/quickguide.md b/docs/src/quickguide.md index f9d07dbc..fadae85c 100644 --- a/docs/src/quickguide.md +++ b/docs/src/quickguide.md @@ -82,7 +82,7 @@ result["objective"] The `"solution"` field contains detailed information about the solution produced by the `run` method. For example, the following dictionary comprehension can be used to inspect the flows in the solution: ``` -Dict(name => data["q"] for (name, data) in result["solution"]["pipe_des"]) +Dict(name => data["q"] for (name, data) in result["solution"]["des_pipe"]) ``` For more information about WaterModels result data see the [WaterModels Result Data Format](@ref) section. @@ -95,13 +95,9 @@ For example, the full mixed-integer nonconvex formulation for design (NLP) can b import KNITRO run_des(data, NLPWaterModel, KNITRO.Optimizer) ``` -and the mixed-integer convex formulation (MICP) can be solved via -```julia -run_des(data, MICPWaterModel, KNITRO.Optimizer) -``` ## Modifying Network Data -The following example demonstrates one way to perform multiple WaterModels solves while modifing network data in Julia. +The following example demonstrates one way to perform multiple WaterModels solves while modifying network data: ```julia run_des(data, MILPRWaterModel, Cbc.Optimizer, ext=Dict(:pipe_breakpoints=>5)) diff --git a/docs/src/specifications.md b/docs/src/specifications.md index c7c47a70..cd80bd81 100644 --- a/docs/src/specifications.md +++ b/docs/src/specifications.md @@ -2,51 +2,49 @@ ## Water Flow (WF) ### Functions -```julia -function_head_loss(wm) -``` ### Objective -```julia -objective_wf(pm) -``` ### Variables -```julia -variable_head(wm, bounded=false) -variable_flow(wm, bounded=false) -variable_volume(wm, bounded=false) -variable_reservoir(wm) -variable_tank(wm) -``` ### Constraints -```julia -for (a, pipe) in ref(wm, :pipe) - constraint_head_loss_pipe(wm, a) -end - -for (a, pump) in ref(wm, :pump) - constraint_head_gain_pump(wm, a, force_on=true) -end - -for (i, node) in ref(wm, :node) - constraint_flow_conservation(wm, i) -end - -for (i, reservoir) in ref(wm, :reservoir) - constraint_source_head(wm, i) - constraint_source_flow(wm, i) -end - -for (i, junction) in ref(wm, :junction) - if junction["demand"] > 0.0 - constraint_sink_flow(wm, i) - end -end - -for i in ids(wm, :tank) - constraint_link_volume(wm, i) - constraint_tank_state(wm, i) -end -``` + +## Multinetwork Water Flow (MN WF) + +### Functions + +### Objective + +### Variables + +### Constraints + +## Optimal Water Flow (OWF) + +### Functions + +### Objective + +### Variables + +### Constraints + +## Multinetwork Optimal Water Flow (MN OWF) + +### Functions + +### Objective + +### Variables + +### Constraints + +## Optimal Design (DES) + +### Functions + +### Objective + +### Variables + +### Constraints From 4beefac52f3a5f6deab07797ebca8c5feb7eb486 Mon Sep 17 00:00:00 2001 From: Byron Tasseff Date: Sat, 22 Aug 2020 13:00:07 -0600 Subject: [PATCH 3/3] Update README. --- README.md | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/README.md b/README.md index 19a89d69..cd584f32 100644 --- a/README.md +++ b/README.md @@ -9,8 +9,8 @@ The code is engineered to decouple problem specifications (e.g., network design, This decoupling enables the definition of a wide variety of optimization formulations and their comparison on common problem specifications. **Core Problem Specifications** -* Water Flow (`wf` and `wf_mn`) - obtain feasible flows for a network -* Optimal Water Flow (`owf` and `owf_mn`) - minimize the cost of network operation +* Water Flow (`wf` and `mn_wf`) - obtain feasible flows for a network +* Optimal Water Flow (`owf` and `mn_owf`) - minimize the cost of network operation * Network Design (`des`) - minimize the cost of network design **Core Network Formulations** @@ -21,6 +21,7 @@ This decoupling enables the definition of a wide variety of optimization formula ## Documentation The package [documentation](https://lanl-ansi.github.io/WaterModels.jl/latest) includes a [quick start guide](https://lanl-ansi.github.io/WaterModels.jl/latest/quickguide). +Be advised that aside from the quick start guide, documentation is under development and may currently be inaccurate. ## Development Community-driven development and enhancement of WaterModels is welcomed and encouraged.