Merge pull request #155 from odow/od/nlp-expr

Update to JuMP v1.15 nonlinear syntax

Project.toml

@@ -5,7 +5,6 @@ repo = "https://github.com/lanl-ansi/WaterModels.jl"
 version = "0.9.3"
 
 [deps]
-HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b"
 InfrastructureModels = "2030c09a-7f63-5d83-885d-db604e0e9cc0"
 Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
 JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
@@ -15,14 +14,16 @@ PolyhedralRelaxations = "2e741578-48fa-11ea-2d62-b52c946f73a0"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [compat]
-HiGHS = "~1"
+HiGHS = "1"
 InfrastructureModels = "~0.6, ~0.7"
 Interpolations = "~0.14"
+Ipopt = "1"
 JSON = "~0.18, ~0.19, ~0.20, ~0.21"
-JuMP = "~0.22, ~0.23, ~1"
+JuMP = "1.15"
+Juniper = "0.9"
 Memento = "~1"
 PolyhedralRelaxations = "~0.3"
-julia = "1.6 - 1"
+julia = "1.6"
 
 [extras]
 HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b"

docs/src/specifications.md

@@ -1,12 +1,6 @@
 # Problem Specifications
 ## Water Flow (WF)
 
-### Functions
-```julia
-# Create head loss functions, if necessary.
-_function_head_loss(wm)
-```
-
 ### Objective
 ```julia
 objective_wf(wm)
@@ -100,12 +94,6 @@ end
 
 ## Multinetwork Water Flow (MN WF)
 
-### Functions
-```julia
-# Create head loss functions, if necessary.
-_function_head_loss(wm)
-```
-
 ### Objective
 ```julia
 objective_wf(wm)

src/core/function.jl

@@ -2,76 +2,23 @@
 # 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
-
-
-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
-    end
-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
-end
-
-
-function _get_alpha_min_1(wm::AbstractWaterModel)
+function _get_alpha(wm::AbstractWaterModel)
     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.")
-    else
-        return alphas[1] - 1.0
     end
+    return alphas[1]
 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))
+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 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/form/crd.jl

@@ -13,8 +13,8 @@ function constraint_pipe_head_loss(
     dhp, dhn = var(wm, n, :dhp_pipe, a), var(wm, n, :dhn_pipe, a)
 
     # Add relaxed constraints for head loss in the positive and negative directions.
-    c_1 = JuMP.@NLconstraint(wm.model, r * head_loss(qp) <= dhp / L)
-    c_2 = JuMP.@NLconstraint(wm.model, r * head_loss(qn) <= dhn / L)
+    c_1 = JuMP.@constraint(wm.model, r * head_loss(wm, qp) <= dhp / L)
+    c_2 = JuMP.@constraint(wm.model, r * head_loss(wm, qn) <= dhn / L)
 
     # Add linear upper bounds on the above convex relaxations.
     rhs_p = r * JuMP.upper_bound(qp)^(exponent - 1.0) * qp
@@ -35,8 +35,8 @@ function constraint_on_off_des_pipe_head_loss(
     dhp, dhn = var(wm, n, :dhp_des_pipe, a), var(wm, n, :dhn_des_pipe, a)
 
     # Add relaxed constraints for head loss in the positive and negative directions.
-    c_1 = JuMP.@NLconstraint(wm.model, r * head_loss(qp) <= dhp / L)
-    c_2 = JuMP.@NLconstraint(wm.model, r * head_loss(qn) <= dhn / L)
+    c_1 = JuMP.@constraint(wm.model, r * head_loss(wm, qp) <= dhp / L)
+    c_2 = JuMP.@constraint(wm.model, r * head_loss(wm, qn) <= dhn / L)
 
     # Add linear upper bounds on the above convex relaxations.
     rhs_p = r * JuMP.upper_bound(qp)^(exponent - 1.0) * qp
@@ -59,7 +59,7 @@ function constraint_on_off_pump_head_gain(wm::AbstractCRDModel, n::Int, a::Int,
     # Define the (relaxed) head gain relationship for the pump.
     head_curve_func_z = _calc_head_curve_function(ref(wm, n, :pump, a), z)
     c_1 = JuMP.@constraint(wm.model, g <= head_curve_func_z(qp))
-    
+
     # Add a constraint that lower-bounds the head gain variable.
     head_curve_func = _calc_head_curve_function(ref(wm, n, :pump, a))
     qp_ub, qp_lb = JuMP.upper_bound(qp), q_min_forward
@@ -87,4 +87,4 @@ function constraint_on_off_pump_power(wm::AbstractCRDModel, n::Int, a::Int, q_mi
         c = JuMP.@constraint(wm.model, power_expr <= P)
         append!(con(wm, n, :on_off_pump_power)[a], [c])
     end
-end
+end

src/form/nc.jl

@@ -190,8 +190,8 @@ function constraint_pipe_head_loss(
     q, h_i, h_j = var(wm, n, :q_pipe, a), var(wm, n, :h, node_fr), var(wm, n, :h, node_to)
 
     # Add nonconvex constraint for the head loss relationship.
-    c_1 = JuMP.@NLconstraint(wm.model, r * head_loss(q) <= (h_i - h_j) / L)
-    c_2 = JuMP.@NLconstraint(wm.model, r * head_loss(q) >= (h_i - h_j) / L)
+    c_1 = JuMP.@constraint(wm.model, r * head_loss(wm, q) <= (h_i - h_j) / L)
+    c_2 = JuMP.@constraint(wm.model, r * head_loss(wm, q) >= (h_i - h_j) / L)
 
     # Append the :pipe_head_loss constraint array.
     append!(con(wm, n, :pipe_head_loss)[a], [c_1, c_2])
@@ -323,8 +323,8 @@ function constraint_on_off_des_pipe_head_loss(
     q, dh = var(wm, n, :q_des_pipe, a), var(wm, n, :dh_des_pipe, a)
 
     # Add nonconvex constraint for the head loss relationship.
-    c_1 = JuMP.@NLconstraint(wm.model, r * head_loss(q) <= dh / L)
-    c_2 = JuMP.@NLconstraint(wm.model, r * head_loss(q) >= dh / L)
+    c_1 = JuMP.@constraint(wm.model, r * head_loss(wm, q) <= dh / L)
+    c_2 = JuMP.@constraint(wm.model, r * head_loss(wm, q) >= dh / L)
 
     # Append the :pipe_head_loss constraint array.
     append!(con(wm, n, :on_off_des_pipe_head_loss)[a], [c_1, c_2])

src/form/ncd.jl

@@ -338,13 +338,13 @@ 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)
+    c_1 = JuMP.@constraint(wm.model, r * head_loss(wm, qp) <= dhp / L)
+    c_2 = JuMP.@constraint(wm.model, r * head_loss(wm, qp) >= 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.@constraint(wm.model, r * head_loss(wm, qn) <= dhn / L)
+    c_4 = JuMP.@constraint(wm.model, r * head_loss(wm, qn) >= 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])
@@ -490,10 +490,10 @@ function constraint_on_off_des_pipe_head_loss(
     dhp, dhn = var(wm, n, :dhp_des_pipe, a), var(wm, n, :dhn_des_pipe, a)
 
     # Add nonconvex constraint for the head loss relationship.
-    c_1 = JuMP.@NLconstraint(wm.model, r * head_loss(qp) <= dhp / L)
-    c_2 = JuMP.@NLconstraint(wm.model, r * head_loss(qp) >= dhp / L)
-    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_1 = JuMP.@constraint(wm.model, r * head_loss(wm, qp) <= dhp / L)
+    c_2 = JuMP.@constraint(wm.model, r * head_loss(wm, qp) >= dhp / L)
+    c_3 = JuMP.@constraint(wm.model, r * head_loss(wm, qn) <= dhn / L)
+    c_4 = JuMP.@constraint(wm.model, r * head_loss(wm, qn) >= dhn / L)
 
     # Append the :pipe_head_loss constraint array.
     append!(con(wm, n, :on_off_des_pipe_head_loss)[a], [c_1, c_2, c_3, c_4])
@@ -1083,4 +1083,4 @@ function constraint_sink_directionality(
     # Add the sink flow direction constraint.
     c = JuMP.@constraint(wm.model, sum_in - sum_out >= 1.0 - out_length)
     con(wm, n, :node_directionality)[i] = c
-end
+end

src/prob/wf.jl

@@ -14,9 +14,6 @@ end
 
 
 function build_wf(wm::AbstractWaterModel)
-    # Create head loss functions, if necessary.
-    _function_head_loss(wm)
-
     # Physical variables.
     variable_head(wm)
     variable_flow(wm)
@@ -103,9 +100,6 @@ end
 
 
 function build_mn_wf(wm::AbstractWaterModel)
-    # Create head loss functions, if necessary.
-    _function_head_loss(wm)
-
     # Get all network IDs in the multinetwork.
     network_ids = sort(collect(nw_ids(wm)))
 
@@ -236,4 +230,4 @@ function build_mn_wf_switching(wm::AbstractWaterModel)
         # Add constraints on the total number of pump switches.
         constraint_on_off_pump_switch(wm, a, network_ids[2:end-1])
     end
-end
+end

test/common.jl

@@ -1,8 +1,7 @@
 function _choose_solver(wm::AbstractWaterModel, nl_solver::JuMP.MOI.OptimizerWithAttributes, mip_solver::JuMP.MOI.OptimizerWithAttributes)
     if isa(wm, AbstractNonlinearModel)
-        f = Juniper.register(head_loss_args(wm)..., autodiff=false)
         return JuMP.optimizer_with_attributes(Juniper.Optimizer,
-            "nl_solver"=>nl_solver, "registered_functions"=>[f], "log_levels"=>[],
+            "nl_solver"=>nl_solver, "log_levels"=>[],
             "allow_almost_solved_integral" => false)
     else
         return mip_solver
@@ -11,4 +10,4 @@ end
 
 function _is_valid_status(status::JuMP.TerminationStatusCode)
     return status in [OPTIMAL, LOCALLY_SOLVED, ALMOST_LOCALLY_SOLVED]
-end
+end

test/function.jl

@@ -1,9 +1,9 @@
 @testset "src/core/function.jl" begin
-    @testset "_get_alpha_min_1 (with non-matching exponents)" begin
+    @testset "_get_alpha (with non-matching exponents)" begin
         data = WaterModels.parse_file("../test/data/epanet/multinetwork/pipe-hw-lps.inp")
         mn_data = WaterModels.make_multinetwork(data)
         wm = instantiate_model(mn_data, NCWaterModel, build_mn_wf)
         WaterModels.ref(wm, 1)[:alpha] = NaN # Change one of the exponents.
-        @test_throws ErrorException WaterModels._get_alpha_min_1(wm)
+        @test_throws ErrorException WaterModels._get_alpha(wm)
     end
-end
+end