Expansion data file and model for expansion pumps

src/WaterModels.jl

@@ -86,6 +86,7 @@ include("form/crd.jl")
 include("form/la.jl")
 include("form/pwlrd.jl")
 include("form/lrd.jl")
+include("form/scip_constraints_variables.jl")
 
 include("prob/wf.jl")
 include("prob/owf.jl")
@@ -98,6 +99,7 @@ include("util/variable_index.jl")
 include("util/bound_problem.jl")
 include("util/pairwise_cuts.jl")
 include("util/obbt.jl")
+include("util/fixing_variables.jl")
 
 # Deprecated functions.
 include("deprecated.jl")

src/core/data.jl

@@ -72,17 +72,19 @@ function _check_connectivity(data::Dict{String,<:Any})
     end
 
     for comp_type in _LINK_COMPONENTS
-        for (i, comp) in data[comp_type]
-            if !(comp["node_fr"] in node_ids)
-                error_message = "From node $(comp["node_fr"]) in "
-                error_message *= "$(replace(comp_type, "_" => " ")) $(i) is not defined."
-                Memento.error(_LOGGER, error_message)
-            end
+        if(haskey(data,comp_type))
+            for (i, comp) in data[comp_type]
+                if !(comp["node_fr"] in node_ids)
+                    error_message = "From node $(comp["node_fr"]) in "
+                    error_message *= "$(replace(comp_type, "_" => " ")) $(i) is not defined."
+                    Memento.error(_LOGGER, error_message)
+                end
 
-            if !(comp["node_to"] in node_ids)
-                error_message = "To node $(comp["node_to"]) in "
-                error_message *= "$(replace(comp_type, "_" => " ")) $(i) is not defined."
-                Memento.error(_LOGGER, error_message)
+                if !(comp["node_to"] in node_ids)
+                    error_message = "To node $(comp["node_to"]) in "
+                    error_message *= "$(replace(comp_type, "_" => " ")) $(i) is not defined."
+                    Memento.error(_LOGGER, error_message)
+                end
             end
         end
     end
@@ -112,17 +114,19 @@ function _check_status(data::Dict{String,<:Any})
     end
 
     for comp_type in _LINK_COMPONENTS
-        for (i, comp) in data[comp_type]
-            if comp["status"] != STATUS_INACTIVE && !(comp["node_fr"] in active_node_ids)
-                warning_message = "Active $(comp_type) $(i) is connected to inactive "
-                warning_message *= "from node $(comp["node_fr"])."
-                Memento.warn(_LOGGER, warning_message)
-            end
+        if(haskey(data,comp_type))
+            for (i, comp) in data[comp_type]
+                if comp["status"] != STATUS_INACTIVE && !(comp["node_fr"] in active_node_ids)
+                    warning_message = "Active $(comp_type) $(i) is connected to inactive "
+                    warning_message *= "from node $(comp["node_fr"])."
+                    Memento.warn(_LOGGER, warning_message)
+                end
 
-            if comp["status"] != STATUS_INACTIVE && !(comp["node_to"] in active_node_ids)
-                warning_message = "Active $(comp_type) $(i) is connected to inactive "
-                warning_message *= "to node $(comp["node_to"])."
-                Memento.warn(_LOGGER, warning_message)
+                if comp["status"] != STATUS_INACTIVE && !(comp["node_to"] in active_node_ids)
+                    warning_message = "Active $(comp_type) $(i) is connected to inactive "
+                    warning_message *= "to node $(comp["node_to"])."
+                    Memento.warn(_LOGGER, warning_message)
+                end
             end
         end
     end
@@ -512,6 +516,14 @@ function _calc_head_max(data::Dict{String, <:Any})
         head_max = max(head_max, node_to["elevation"] + head_gain)
     end
 
+    for (i, pump) in wm_data["ne_pump"]
+        # Consider possible pump head gains in computation of head_max.
+        node_fr = wm_data["node"][string(pump["node_fr"])]
+        node_to = wm_data["node"][string(pump["node_to"])]
+        head_gain = _calc_pump_head_gain_max(pump, node_fr, node_to)
+        head_max = max(head_max, node_to["elevation"] + head_gain)
+    end
+
     for (i, regulator) in wm_data["regulator"]
         # Consider possible downstream regulator heads in computation of head_max.
         p_setting, node_to_index = regulator["setting"], string(regulator["node_to"])
@@ -1023,7 +1035,7 @@ end
 
 function _apply_pipe_unit_transform!(data::Dict{String,<:Any}, transform_length::Function, head_loss::String)
     wm_data = get_wm_data(data)
-    
+
     if !haskey(wm_data, "pipe")
         return
     end
@@ -1139,6 +1151,78 @@ function _apply_pump_unit_transform!(
     end
 end
 
+function _apply_ne_pump_unit_transform!(
+    data::Dict{String,<:Any}, transform_mass::Function, transform_flow::Function,
+    transform_length::Function, transform_time::Function)
+    wm_data = get_wm_data(data)
+
+    if !haskey(wm_data, "ne_pump")
+        return
+    end
+
+    power_scalar = transform_mass(1.0) * transform_length(1.0)^2 / transform_time(1.0)^3
+    energy_scalar = transform_mass(1.0) * transform_length(1.0)^2 / transform_time(1.0)^2
+
+    for (_, ne_pump) in wm_data["ne_pump"]
+        if haskey(ne_pump, "head_curve")
+            ne_pump["head_curve"] = [(transform_flow(x[1]), x[2]) for x in ne_pump["head_curve"]]
+            ne_pump["head_curve"] = [(x[1], transform_length(x[2])) for x in ne_pump["head_curve"]]
+        end
+
+        if haskey(ne_pump, "efficiency_curve")
+            ne_pump["efficiency_curve"] = [(transform_flow(x[1]), x[2]) for x in ne_pump["efficiency_curve"]]
+        end
+
+        if haskey(ne_pump, "energy_price")
+            ne_pump["energy_price"] /= energy_scalar
+        end
+
+        if haskey(ne_pump, "E")
+            ne_pump["E"] *= energy_scalar
+        end
+
+        if haskey(ne_pump, "power_fixed")
+            ne_pump["power_fixed"] *= power_scalar
+        end
+
+        if haskey(ne_pump, "P")
+            ne_pump["P"] *= power_scalar
+        end
+
+        if haskey(ne_pump, "c")
+            ne_pump["c"] *= energy_scalar
+        end
+
+        if haskey(ne_pump, "min_inactive_time")
+            ne_pump["min_inactive_time"] = transform_time(ne_pump["min_inactive_time"])
+        end
+
+        if haskey(ne_pump, "min_active_time")
+            ne_pump["min_active_time"] = transform_time(ne_pump["min_active_time"])
+        end
+
+        if haskey(ne_pump, "power_per_unit_flow")
+            ne_pump["power_per_unit_flow"] *= power_scalar / transform_flow(1.0)
+        end
+    end
+
+    if haskey(wm_data, "time_series") && haskey(wm_data["time_series"], "ne_pump")
+        for ne_pump in values(wm_data["time_series"]["ne_pump"])
+            if haskey(ne_pump, "energy_price")
+                ne_pump["energy_price"] ./= energy_scalar
+            end
+
+            if haskey(ne_pump, "power_fixed")
+                ne_pump["power_fixed"] .*= power_scalar
+            end
+
+            if haskey(ne_pump, "power_per_unit_flow")
+                ne_pump["power_per_unit_flow"] .*= power_scalar / transform_flow(1.0)
+            end
+        end
+    end
+end
+
 
 function _apply_regulator_unit_transform!(data::Dict{String,<:Any}, transform_head::Function, transform_length::Function)
     wm_data = get_wm_data(data)
@@ -1169,7 +1253,7 @@ function _calc_scaled_gravity(data::Dict{String, <:Any})
 
     if wm_data["per_unit"]
         base_time = 1.0 / _calc_time_per_unit_transform(wm_data)(1.0)
-        base_length = 1.0 / _calc_length_per_unit_transform(wm_data)(1.0) 
+        base_length = 1.0 / _calc_length_per_unit_transform(wm_data)(1.0)
         return _calc_scaled_gravity(base_length, base_time)
     else
         return _GRAVITY
@@ -1368,6 +1452,8 @@ function _transform_nw!(
     _apply_des_pipe_unit_transform!(wm_nw_data, length_transform, head_loss)
     _apply_pump_unit_transform!(wm_nw_data, mass_transform,
         flow_transform, length_transform, time_transform)
+    _apply_ne_pump_unit_transform!(wm_nw_data, mass_transform,
+        flow_transform, length_transform, time_transform)
     _apply_regulator_unit_transform!(wm_nw_data, head_transform, length_transform)
 
     # Apply transformations to nodal components.

src/core/function.jl

@@ -1,84 +1,27 @@
-#################################################################################
-# This file defines the nonlinear head loss functions for water systems models. #
-#################################################################################
 
 
-function _f_alpha(alpha::Float64; convex::Bool = false)
-    if convex
-        return function (x::Float64)
-            return (x * x)^(0.5 + 0.5 * alpha)
-        end
-    else
-        return function (x::Float64)
-            return sign(x) * (x * x)^(0.5 + 0.5 * alpha)
-        end
-    end
-end
-
+# This file defines the nonlinear head loss functions for water systems models. #
+#################################################################################
 
-function _df_alpha(alpha::Float64; convex::Bool = false)
-    if convex
-        return function (x::Float64)
-            return (1.0 + alpha) * sign(x) * (x * x)^(0.5 * alpha)
-        end
-    else
-        return function (x::Float64)
-            return (1.0 + alpha) * (x * x)^(0.5 * alpha)
-        end
+function _get_alpha(wm::AbstractWaterModel)
+    head_loss_form = wm.ref[:it][wm_it_sym][:head_loss]
+    alphas = uppercase(head_loss_form) == "H-W" ? 1.852 : 2.0
+    # alphas = [ref(wm, nw, :alpha) for nw in nw_ids(wm)]
+    if any(x -> x != alphas[1], alphas)
+        Memento.error(_LOGGER, "Head loss exponents are different across the multinetwork.")
     end
+    return alphas[1]
 end
 
-
-function _d2f_alpha(alpha::Float64; convex::Bool = false)
-    if convex
-        return function (x::Float64)
-            return alpha * (1.0 + alpha) * (x * x)^(0.5 * alpha - 0.5)
-        end
-    else
-        return function (x::Float64)
-            return x != 0.0 ?
-                   sign(x) * alpha * (1.0 + alpha) * (x * x)^(0.5 * alpha - 0.5) :
-                   prevfloat(Inf)
-        end
-    end
+function head_loss(wm::Union{NCDWaterModel,CRDWaterModel}, x)
+    p = _get_alpha(wm)
+    # An expression equivalent to abspower(x, p) = abs(x)^p
+    # return JuMP.@expression(wm.model, ifelse(x > 0, x^p, (-x)^p))
+    return JuMP.@expression(wm.model, (x^2)^(p / 2))
 end
 
-
-function _get_alpha_min_1(wm::AbstractWaterModel)
-    return _get_exponent_from_head_loss_form(wm.ref[:it][wm_it_sym][:head_loss]) - 1.0
-end
-
-
-function head_loss_args(wm::Union{NCDWaterModel,CRDWaterModel})
-    alpha_m1 = _get_alpha_min_1(wm)
-    return (
-        :head_loss,
-        1,
-        _f_alpha(alpha_m1, convex = true),
-        _df_alpha(alpha_m1, convex = true),
-        _d2f_alpha(alpha_m1, convex = true),
-    )
-end
-
-
-function head_loss_args(wm::NCWaterModel)
-    alpha_m1 = _get_alpha_min_1(wm)
-
-    return (
-        :head_loss,
-        1,
-        _f_alpha(alpha_m1, convex = false),
-        _df_alpha(alpha_m1, convex = false),
-        _d2f_alpha(alpha_m1, convex = false),
-    )
-end
-
-
-function _function_head_loss(wm::AbstractWaterModel)
-    # By default, head loss is not defined by nonlinear registered functions.
-end
-
-
-function _function_head_loss(wm::AbstractNonlinearModel)
-    JuMP.register(wm.model, head_loss_args(wm)...)
+function head_loss(wm::AbstractNonlinearModel, x)
+    p = _get_alpha(wm)
+    # An expression equivalent to signpower(x, p) = sign(x) * abs(x)^p
+    return JuMP.@expression(wm.model, ifelse(x > 0, x^p, -(-x)^p))
 end

src/core/objective.jl

@@ -85,32 +85,25 @@ cost of building and operating all network expansion components is minimized.
 """
 function objective_ne(wm::AbstractWaterModel)::JuMP.AffExpr
     # Get all network IDs in the multinetwork.
-    network_ids = sort(collect(nw_ids(wm)))
+    first_network_id = sort(collect(nw_ids(wm)))[1]
 
-    # Find the network IDs over which the objective will be defined.
-    if length(network_ids) > 1
-        network_ids_flow = network_ids[1:end-1]
-    else
-        network_ids_flow = network_ids
-    end
 
     # Initialize the objective expression to zero.
     objective = JuMP.AffExpr(0.0)
 
-    for n in network_ids_flow
-        # Get the set of network expansion short pipes at time index `n`.
-        for (a, ne_short_pipe) in ref(wm, n, :ne_short_pipe)
-            # Add the cost of network expansion component `a` at time period `n`.
-            term = ne_short_pipe["construction_cost"] * var(wm, n, :z_ne_short_pipe, a)
-            JuMP.add_to_expression!(objective, term)
-        end
+n = first_network_id
+    # Get the set of network expansion short pipes at time index `n`.
+    for (a, ne_short_pipe) in ref(wm, n, :ne_short_pipe)
+        # Add the cost of network expansion component `a` at time period `n`.
+        term = ne_short_pipe["construction_cost"] * var(wm, n, :z_ne_short_pipe, a)
+        JuMP.add_to_expression!(objective, term)
+    end
 
-        # Get the set of network expansion pumps at time index `n`.
-        for (a, ne_pump) in ref(wm, n, :ne_pump)
-            # Add the cost of network expansion component `a` at time period `n`.
-            term = ne_pump["construction_cost"] * var(wm, n, :z_ne_pump, a)
-            JuMP.add_to_expression!(objective, term)
-        end
+    # Get the set of network expansion pumps at time index `n`.
+    for (a, ne_pump) in ref(wm, n, :ne_pump)
+        # Add the cost of network expansion component `a` at time period `n`.
+        term = ne_pump["construction_cost"] * var(wm, n, :x_ne_pump, a)
+        JuMP.add_to_expression!(objective, term)
     end
 
     # Minimize the total cost of network expansion.
@@ -180,3 +173,42 @@ function objective_owf(wm::AbstractWaterModel)::JuMP.AffExpr
     # Minimize the cost of network operation.
     return JuMP.@objective(wm.model, JuMP.MIN_SENSE, objective)
 end
+
+
+
+"""
+    objective_max_demand(wm::AbstractWaterModel)::JuMP.AffExpr
+"""
+function objective_max_scaled_demand(wm::AbstractWaterModel)::JuMP.AffExpr
+    # Get all network IDs in the multinetwork.
+    network_ids = sort(collect(nw_ids(wm)))
+
+    # Find the network IDs over which the objective will be defined.
+    if length(network_ids) > 1
+        network_ids_flow = network_ids[1:end-1]
+    else
+        network_ids_flow = network_ids
+    end
+
+    # Initialize the objective expression to zero.
+    objective = JuMP.AffExpr(0.0)
+
+    scale = 0.0
+    for n in network_ids_flow
+        # Get the set of dispatchable demands at time index `n`.
+        dispatchable_demands = ref(wm, n, :dispatchable_demand)
+
+        for (i, demand) in filter(x -> x.second["flow_min"] >= 0.0, dispatchable_demands)
+            # Add the volume delivered at demand `i` and time period `n`.
+            # coeff = get(demand, "priority", 1.0) * ref(wm, n, :time_step)
+            scale = max(scale,demand["flow_max"])
+            JuMP.add_to_expression!(objective, var(wm, n, :q_demand, i))
+        end
+    end
+    if(scale==0.0)
+        scale = 1.0
+    end
+    println("Using this scale for water objective $(scale)*****")
+    # Maximize the total amount of prioritized water volume delivered.
+    return JuMP.@objective(wm.model, JuMP.MAX_SENSE, objective/scale)
+end