head_loss_function

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

@@ -173,3 +173,39 @@ 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 = 1.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
+    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

src/form/ncd.jl

@@ -343,13 +343,15 @@ function constraint_pipe_head_loss(
 )
     # Add constraints for positive flow and head difference.
     qp, dhp = var(wm, n, :qp_pipe, a), var(wm, n, :dhp_pipe, a)
-    c_1 = JuMP.@NLconstraint(wm.model, r * head_loss(qp) <= dhp / L)
-    c_2 = JuMP.@NLconstraint(wm.model, r * head_loss(qp) >= dhp / L)
+    head_loss_form = wm.ref[:it][wm_it_sym][:head_loss]
+    p = uppercase(head_loss_form) == "H-W" ? 1.852 : 2.0
+    c_1 = JuMP.@NLconstraint(wm.model, r * (qp^p) <= dhp / L)
+    c_2 = JuMP.@NLconstraint(wm.model, r * (qp^p) >= dhp / L)
 
     # Add constraints for negative flow and head difference.
     qn, dhn = var(wm, n, :qn_pipe, a), var(wm, n, :dhn_pipe, a)
-    c_3 = JuMP.@NLconstraint(wm.model, r * head_loss(qn) <= dhn / L)
-    c_4 = JuMP.@NLconstraint(wm.model, r * head_loss(qn) >= dhn / L)
+    c_3 = JuMP.@NLconstraint(wm.model, r * (qn^p) <= dhn / L)
+    c_4 = JuMP.@NLconstraint(wm.model, r * (qn^p) >= dhn / L)
 
     # Append the :pipe_head_loss constraint array.
     append!(con(wm, n, :pipe_head_loss)[a], [c_1, c_2, c_3, c_4])

src/form/scip_constraints_variables.jl

@@ -59,112 +59,119 @@ function constraint_pipe_head_loss_scip_lp(
     q_max_reverse::Float64,
     q_min_forward::Float64,
 )
-    # Get the variable for flow directionality.
-    y = var(wm, n, :y_pipe, a)
-
-    # Get variables for positive flow and head difference.
-    qp = var(wm, n, :qp_pipe, a)
-    dhp = var(wm, n, :dhp_pipe, a)
-
-    # Get the corresponding positive flow partitioning.
-    partition_p = get_pipe_flow_partition_positive(ref(wm, n, :pipe, a))
-
-    # Loop over consequential points (i.e., those that have nonzero head loss).
-    for flow_value in filter(x -> x > 0.0, partition_p)
-        # Add a linear outer approximation of the convex relaxation at `flow_value`.
-        lhs = r * _calc_head_loss_oa(qp, y, flow_value, exponent)
-
-        if minimum(abs.(lhs.terms.vals)) >= 1.0e-4
-            # Compute a scaling factor to normalize the constraint.
-            scalar = _get_scaling_factor(vcat(lhs.terms.vals, [1.0 / L]))
-
-            # Add outer-approximation of the head loss constraint.
-            c = JuMP.@constraint(wm.model, scalar * lhs <= scalar * dhp / L)
-
-            # Append the :pipe_head_loss constraint array.
-            append!(con(wm, n, :pipe_head_loss)[a], [c])
-        end
-    end
-
-    # Get the corresponding min/max positive directed flows (when active).
-    qp_min_forward, qp_max = max(0.0, q_min_forward), maximum(partition_p)
-
-    if qp_min_forward != qp_max
-        # Compute scaled version of the head loss overapproximation.
-        f_1, f_2 = r * qp_min_forward^exponent, r * qp_max^exponent
-        f_slope = (f_2 - f_1) / (qp_max - qp_min_forward)
-        f_lb_line = f_slope * (qp - qp_min_forward * y) + f_1 * y
-
-        # Compute a scaling factor to normalize the constraint.
-        scalar = _get_scaling_factor(vcat(f_lb_line.terms.vals, [1.0 / L]))
-
-        # Add upper-bounding line of the head loss constraint.
-        c = JuMP.@constraint(wm.model, scalar * dhp / L <= scalar * f_lb_line)
-
-        # Append the :pipe_head_loss constraint array.
-        append!(con(wm, n, :pipe_head_loss)[a], [c])
-    elseif qp_max == 0.0
-        c = JuMP.@constraint(wm.model, dhp == 0.0)
-        append!(con(wm, n, :pipe_head_loss)[a], [c])
-    else
-        f_q = r * qp_max^exponent
-        scalar = _get_scaling_factor([f_q == 0.0 ? 1.0 : f_q, 1.0 / L])
-        c = JuMP.@constraint(wm.model, scalar * dhp / L == scalar * f_q * y)
-        append!(con(wm, n, :pipe_head_loss)[a], [c])
-    end
-
-    # Get variables for negative flow and head difference.
-    qn = var(wm, n, :qn_pipe, a)
-    dhn = var(wm, n, :dhn_pipe, a)
-
-    # Get the corresponding negative flow partitioning (negated).
-    partition_n = sort(-get_pipe_flow_partition_negative(ref(wm, n, :pipe, a)))
-
-    # Loop over consequential points (i.e., those that have nonzero head loss).
-    for flow_value in filter(x -> x > 0.0, partition_n)
-        # Add a linear outer approximation of the convex relaxation at `flow_value`.
-        lhs = r * _calc_head_loss_oa(qn, 1.0 - y, flow_value, exponent)
+# Get the variable for flow directionality.
+y = var(wm, n, :y_pipe, a)
+
+# Get variables for positive flow and head difference.
+qp, dhp = var(wm, n, :qp_pipe, a), var(wm, n, :dhp_pipe, a)
+
+# Get positively-directed convex combination variables.
+lambda_p, x_p = var(wm, n, :lambda_p_pipe), var(wm, n, :x_p_pipe)
+
+# Get the corresponding positive flow partitioning.
+partition_p = get_pipe_flow_partition_positive(ref(wm, n, :pipe, a))
+bp_range, bp_range_m1 = 1:length(partition_p), 1:length(partition_p)-1
+
+# Add constraints for the positive flow piecewise approximation.
+c_1 = JuMP.@constraint(wm.model, sum(lambda_p[a, k] for k in bp_range) == y)
+qp_sum = sum(partition_p[k] * lambda_p[a, k] for k in bp_range)
+scalar = _get_scaling_factor(vcat(qp_sum.terms.vals, [1.0]))
+c_2 = JuMP.@constraint(wm.model, scalar * qp_sum == scalar * qp)
+append!(con(wm, n, :pipe_head_loss, a), [c_1, c_2])
+
+if length(partition_p) > 1
+    # If there are multiple points, constrain the convex combination.
+    c_3 = JuMP.@constraint(wm.model, sum(x_p[a, k] for k in bp_range_m1) == y)
+    c_4 = JuMP.@constraint(wm.model, lambda_p[a, 1] <= x_p[a, 1])
+    c_5 = JuMP.@constraint(wm.model, lambda_p[a, bp_range[end]] <= x_p[a, bp_range_m1[end]])
+    append!(con(wm, n, :pipe_head_loss, a), [c_3, c_4, c_5])
+end
 
-        if minimum(abs.(lhs.terms.vals)) >= 1.0e-4
-            # Compute a scaling factor to normalize the constraint.
-            scalar = _get_scaling_factor(vcat(lhs.terms.vals, [1.0 / L]))
+# Add a constraint that upper-bounds the head loss variable.
+if maximum(partition_p) != 0.0
+    f_p = r * partition_p.^exponent
+    f_p_ub_expr = sum(f_p[k] * lambda_p[a, k] for k in bp_range)
+    scalar = _get_scaling_factor(vcat(f_p_ub_expr.terms.vals, [1.0 / L]))
+    c_6 = JuMP.@constraint(wm.model, scalar * dhp / L <= scalar * f_p_ub_expr)
+    append!(con(wm, n, :pipe_head_loss, a), [c_6])
+else
+    c_6 = JuMP.@constraint(wm.model, dhp == 0.0)
+    append!(con(wm, n, :pipe_head_loss, a), [c_6])
+end
 
-            # Add outer-approximation of the head loss constraint.
-            c = JuMP.@constraint(wm.model, scalar * lhs <= scalar * dhn / L)
+# Loop over consequential points (i.e., those that have nonzero head loss).
+for flow_value in filter(x -> x > 0.0, partition_p)
+    # Add head loss outer (i.e., lower) approximations.
+    lhs_p = r * _calc_head_loss_oa(qp, y, flow_value, exponent)
 
-            # Append the :pipe_head_loss constraint array.
-            append!(con(wm, n, :pipe_head_loss)[a], [c])
-        end
+    if minimum(abs.(lhs_p.terms.vals)) >= 1.0e-4
+        scalar = _get_scaling_factor(vcat(lhs_p.terms.vals, [1.0 / L]))
+        c_7 = JuMP.@constraint(wm.model, scalar * lhs_p <= scalar * dhp / L)
+        append!(con(wm, n, :pipe_head_loss, a), [c_7])
     end
+end
 
-    # Get the corresponding maximum negative directed flow (when active).
-    qn_min_forward, qn_max = max(0.0, -q_max_reverse), maximum(partition_n)
+for k in 2:length(partition_p)-1
+    # Add the adjacency constraints for piecewise variables.
+    c_8 = JuMP.@constraint(wm.model, lambda_p[a, k] <= x_p[a, k-1] + x_p[a, k])
+    append!(con(wm, n, :pipe_head_loss, a), [c_8])
+end
 
-    if qn_min_forward != qn_max
-        # Compute scaled version of the head loss overapproximation.
-        f_1, f_2 = r * qn_min_forward^exponent, r * qn_max^exponent
-        f_slope = (f_2 - f_1) / (qn_max - qn_min_forward)
-        f_lb_line = f_slope * (qn - qn_min_forward * (1.0 - y)) + f_1 * (1.0 - y)
+# Get variables for negative flow and head difference.
+qn, dhn = var(wm, n, :qn_pipe, a), var(wm, n, :dhn_pipe, a)
+
+# Get negatively-directed convex combination variable.
+lambda_n, x_n = var(wm, n, :lambda_n_pipe), var(wm, n, :x_n_pipe)
+
+# Get the corresponding negative flow partitioning (negated).
+partition_n = sort(-get_pipe_flow_partition_negative(ref(wm, n, :pipe, a)))
+bn_range, bn_range_m1 = 1:length(partition_n), 1:length(partition_n)-1
+
+# Add constraints for the negative flow piecewise approximation.
+c_9 = JuMP.@constraint(wm.model, sum(lambda_n[a, k] for k in bn_range) == 1.0 - y)
+qn_sum = sum(partition_n[k] * lambda_n[a, k] for k in bn_range)
+scalar = _get_scaling_factor(vcat(qn_sum.terms.vals, [1.0]))
+c_10 = JuMP.@constraint(wm.model, scalar * qn_sum == scalar * qn)
+append!(con(wm, n, :pipe_head_loss, a), [c_9, c_10])
+
+if length(partition_n) > 1
+    # If there are multiple points, constrain the convex combination.
+    c_11 = JuMP.@constraint(wm.model, sum(x_n[a, k] for k in bn_range_m1) == 1.0 - y)
+    c_12 = JuMP.@constraint(wm.model, lambda_n[a, 1] <= x_n[a, 1])
+    c_13 = JuMP.@constraint(wm.model, lambda_n[a, bn_range[end]] <= x_n[a, bn_range_m1[end]])
+    append!(con(wm, n, :pipe_head_loss, a), [c_11, c_12, c_13])
+end
 
-        # Compute a scaling factor to normalize the constraint.
-        scalar = _get_scaling_factor(vcat(f_lb_line.terms.vals, [1.0 / L]))
+# Add a constraint that upper-bounds the head loss variable.
+if maximum(partition_n) != 0.0
+    f_n = r .* partition_n.^exponent
+    f_n_ub_expr = sum(f_n[k] * lambda_n[a, k] for k in bn_range)
+    scalar = _get_scaling_factor(vcat(f_n_ub_expr.terms.vals, [1.0 / L]))
+    c_14 = JuMP.@constraint(wm.model, scalar * dhn / L <= scalar * f_n_ub_expr)
+    append!(con(wm, n, :pipe_head_loss, a), [c_14])
+else
+    c_14 = JuMP.@constraint(wm.model, dhn == 0.0)
+    append!(con(wm, n, :pipe_head_loss, a), [c_14])
+end
 
-        # Add upper-bounding line of the head loss constraint.
-        c = JuMP.@constraint(wm.model, scalar * dhn / L <= scalar * f_lb_line)
+# Loop over consequential points (i.e., those that have nonzero head loss).
+for flow_value in filter(x -> x > 0.0, partition_n)
+    # Add head loss outer (i.e., lower) approximations.
+    lhs_n = r * _calc_head_loss_oa(qn, 1.0 - y, flow_value, exponent)
 
-        # Append the :pipe_head_loss constraint array.
-        append!(con(wm, n, :pipe_head_loss)[a], [c])
-    elseif qn_max == 0.0
-        c = JuMP.@constraint(wm.model, dhn == 0.0)
-        append!(con(wm, n, :pipe_head_loss)[a], [c])
-    else
-        f_q = r * qn_max^exponent
-        scalar = _get_scaling_factor([f_q == 0.0 ? 1.0 : f_q, 1.0 / L])
-        c = JuMP.@constraint(wm.model, scalar * dhn / L == scalar * f_q * (1.0 - y))
-        append!(con(wm, n, :pipe_head_loss)[a], [c])
+    if minimum(abs.(lhs_n.terms.vals)) >= 1.0e-4
+        scalar = _get_scaling_factor(vcat(lhs_n.terms.vals, [1.0 / L]))
+        c_15 = JuMP.@constraint(wm.model, scalar * lhs_n <= scalar * dhn / L)
+        append!(con(wm, n, :pipe_head_loss, a), [c_15])
     end
 end
 
+for k in 2:length(partition_n)-1
+    # Add the adjacency constraints for piecewise variables.
+    c_16 = JuMP.@constraint(wm.model, lambda_n[a, k] <= x_n[a, k-1] + x_n[a, k])
+    append!(con(wm, n, :pipe_head_loss, a), [c_16])
+end
+end