Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Update to JuMP v1.15 nonlinear syntax #155

Merged
merged 7 commits into from
Sep 17, 2023
Merged
Show file tree
Hide file tree
Changes from 5 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
10 changes: 6 additions & 4 deletions Project.toml
Original file line number Diff line number Diff line change
Expand Up @@ -5,24 +5,26 @@ 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"
JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
Memento = "f28f55f0-a522-5efc-85c2-fe41dfb9b2d9"
Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
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"
Juniper = "0.9"
Memento = "~1"
PolyhedralRelaxations = "~0.3"
julia = "1.6 - 1"
julia = "1.6"

[extras]
HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b"
Expand Down
12 changes: 0 additions & 12 deletions docs/src/specifications.md
Original file line number Diff line number Diff line change
@@ -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)
Expand Down Expand Up @@ -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)
Expand Down
75 changes: 11 additions & 64 deletions src/core/function.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
12 changes: 6 additions & 6 deletions src/form/crd.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
odow marked this conversation as resolved.
Show resolved Hide resolved
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
Expand All @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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
8 changes: 4 additions & 4 deletions src/form/nc.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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])
Expand Down Expand Up @@ -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])
Expand Down
18 changes: 9 additions & 9 deletions src/form/ncd.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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])
Expand Down Expand Up @@ -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])
Expand Down Expand Up @@ -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
8 changes: 1 addition & 7 deletions src/prob/wf.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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)))

Expand Down Expand Up @@ -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
5 changes: 2 additions & 3 deletions test/common.jl
Original file line number Diff line number Diff line change
@@ -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"=>[],
odow marked this conversation as resolved.
Show resolved Hide resolved
"allow_almost_solved_integral" => false)
else
return mip_solver
Expand All @@ -11,4 +10,4 @@ end

function _is_valid_status(status::JuMP.TerminationStatusCode)
return status in [OPTIMAL, LOCALLY_SOLVED, ALMOST_LOCALLY_SOLVED]
end
end
6 changes: 3 additions & 3 deletions test/function.jl
Original file line number Diff line number Diff line change
@@ -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
3 changes: 3 additions & 0 deletions test/runtests.jl
Original file line number Diff line number Diff line change
@@ -1,3 +1,6 @@
import Pkg
Pkg.pkg"add JuMP#master Statistics"
odow marked this conversation as resolved.
Show resolved Hide resolved

using WaterModels

import HiGHS
Expand Down