Skip to content

Commit

Permalink
initial updates for nl interface (#256)
Browse files Browse the repository at this point in the history
  • Loading branch information
ccoffrin authored Jan 31, 2024
1 parent 70ae501 commit 7a0dc25
Show file tree
Hide file tree
Showing 8 changed files with 34 additions and 33 deletions.
6 changes: 5 additions & 1 deletion CHANGELOG.md
Original file line number Diff line number Diff line change
@@ -1,7 +1,11 @@
# GasModels.jl Change Log

## Staged
## Staged
- nothing

## v0.10.0
- update to new JuMP nonlinear interface (breaking)
- drop support for JuMP v0.22 and v0.23 (breaking)

## v0.9.3
- update Memento dependency version
Expand Down
13 changes: 5 additions & 8 deletions Project.toml
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@ name = "GasModels"
uuid = "5e113713-6c35-5477-b766-e1109486666f"
authors = ["Russell Bent <[email protected]>", "Kaarthik Sundar <[email protected]>", "David Fobes <[email protected]>"]
repo = "https://github.com/lanl-ansi/GasModels.jl.git"
version = "0.9.3"
version = "0.10.0"

[deps]
Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
Expand All @@ -11,35 +11,32 @@ InfrastructureModels = "2030c09a-7f63-5d83-885d-db604e0e9cc0"
JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
Memento = "f28f55f0-a522-5efc-85c2-fe41dfb9b2d9"
PolyhedralRelaxations = "2e741578-48fa-11ea-2d62-b52c946f73a0"
Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
XMLDict = "228000da-037f-5747-90a9-8195ccbf91a5"
ZipFile = "a5390f91-8eb1-5f08-bee0-b1d1ffed6cea"
PolyhedralRelaxations = "2e741578-48fa-11ea-2d62-b52c946f73a0"



[compat]
HiGHS = "~0.3, ~1"
Dierckx = "~0.4, ~0.5"
HiGHS = "~0.3, ~1"
InfrastructureModels = "~0.6, ~0.7"
Ipopt = "~0.8, ~0.9, ~1"
JSON = "~0.18, ~0.19, ~0.20, ~0.21"
JuMP = "~0.22, ~0.23, 1"
JuMP = "1.15"
Juniper = ">= 0.4"
Memento = "~1.0, ~1.1, ~1.2, ~1.3, ~1.4"
PolyhedralRelaxations = "~0.3"
XMLDict = "~0.4"
ZipFile = "~0.8, ~0.9, ~0.10"
julia = "^1"


[extras]
HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b"
MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9"
Juniper = "2ddba703-00a4-53a7-87a5-e8b9971dde84"
Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[targets]
Expand Down
8 changes: 4 additions & 4 deletions src/core/constraint_transient.jl
Original file line number Diff line number Diff line change
Expand Up @@ -30,7 +30,7 @@ function constraint_pipe_momentum_balance(gm::AbstractGasModel, nw::Int, pipe_id
p_fr = var(gm, nw, :density, fr_junction)
p_to = var(gm, nw, :density, to_junction)
f = var(gm, nw, :pipe_flux_avg, pipe_id)
_add_constraint!(gm, nw, :pipe_momentum_balance, pipe_id, JuMP.@NLconstraint(gm.model, p_fr^2 - p_to^2 - resistance * f * abs(f) == 0))
_add_constraint!(gm, nw, :pipe_momentum_balance, pipe_id, JuMP.@constraint(gm.model, p_fr^2 - p_to^2 - resistance * f * abs(f) == 0))
end

"Constraint: inclined pipe momentum balance.
Expand All @@ -46,7 +46,7 @@ function constraint_inclined_pipe_momentum_balance(gm::AbstractGasModel, nw::Int

r_1 = resistance_1
r_2 = resistance_2
_add_constraint!(gm, nw, :inclined_pipe_momentum_balance, pipe_id, JuMP.@NLconstraint(gm.model, p_to^2 - exp(r_2)*p_fr^2 == r_1*(exp(r_2) - 1)*f*abs(f)))
_add_constraint!(gm, nw, :inclined_pipe_momentum_balance, pipe_id, JuMP.@constraint(gm.model, p_to^2 - exp(r_2)*p_fr^2 == r_1*(exp(r_2) - 1)*f*abs(f)))
end

"Constraint: non-slack junction mass balance"
Expand All @@ -59,7 +59,7 @@ function constraint_pipe_physics_ideal(gm::AbstractGasModel, nw::Int, pipe_id::I
p_fr = var(gm, nw, :density, fr_junction)
p_to = var(gm, nw, :density, to_junction)
f = var(gm, nw, :pipe_flux, pipe_id)
_add_constraint!(gm, nw, :pipe_physics_ideal, pipe_id, JuMP.@NLconstraint(gm.model, p_fr^2 - p_to^2 - resistance * f * abs(f) == 0))
_add_constraint!(gm, nw, :pipe_physics_ideal, pipe_id, JuMP.@constraint(gm.model, p_fr^2 - p_to^2 - resistance * f * abs(f) == 0))
end

"Constraint: aggregate withdrawal at transfer points computation"
Expand Down Expand Up @@ -110,7 +110,7 @@ end

"Constraint: compressor power"
function constraint_compressor_power(gm::AbstractGasModel, nw::Int, compressor_id::Int, compressor_power_expr, compressor_power_var)
_add_constraint!(gm, nw, :compressor_power, compressor_id, JuMP.@NLconstraint(gm.model, compressor_power_var == compressor_power_expr))
_add_constraint!(gm, nw, :compressor_power, compressor_id, JuMP.@constraint(gm.model, compressor_power_var == compressor_power_expr))
end


Expand Down
6 changes: 3 additions & 3 deletions src/core/objective.jl
Original file line number Diff line number Diff line change
Expand Up @@ -64,8 +64,8 @@ function objective_min_compressor_energy(gm::AbstractGasModel, nws = [nw_id_defa

# some solvers only support support nonlinear objectives by placing them in the constraints
z = JuMP.@variable(gm.model)
JuMP.@NLconstraint(gm.model, z >= sum(sum((r[n][i]^m - 1) * f[n][i] for (i, compressor) in ref(gm, n, :compressor)) for n in nws))
return JuMP.@NLobjective(gm.model, Min, z)
JuMP.@constraint(gm.model, z >= sum(sum((r[n][i]^m - 1) * f[n][i] for (i, compressor) in ref(gm, n, :compressor)) for n in nws))
return JuMP.@objective(gm.model, Min, z)
end


Expand Down Expand Up @@ -119,7 +119,7 @@ function objective_min_economic_costs(gm::AbstractGasModel, nws = [nw_id_default

# prices are already normalized by base_flow, so we also need to normalize compressor power by base_flow in the objective
z = JuMP.@variable(gm.model)
JuMP.@NLconstraint(gm.model, z >= sum(
JuMP.@constraint(gm.model, z >= sum(
economic_weighting * sum(-load_prices[n][i] * fl[n][i] for i in load_set[n]) +
economic_weighting *
sum(-transfer_prices[n][i] * ft[n][i] for i in transfer_set[n]) +
Expand Down
2 changes: 1 addition & 1 deletion src/core/storage_archived.jl
Original file line number Diff line number Diff line change
Expand Up @@ -366,7 +366,7 @@ function constraint_storage_well_momentum_balance(gm::AbstractGasModel, nw::Int,
rho_top = var(gm, nw, :well_density, storage_id)[i]
rho_bottom = var(gm, nw, :well_density, storage_id)[i+1]
phi_avg = var(gm, nw, :well_flux_avg, storage_id)[i]
GasModels._add_constraint!(gm, nw, :well_ideal_momentum_balance, storage_id * 1000 + i, JuMP.@NLconstraint(gm.model, exp(beta) * rho_bottom^2 - rho_top^2 == (-resistance * phi_avg * abs(phi_avg)) * (exp(beta) - 1) / beta))
GasModels._add_constraint!(gm, nw, :well_ideal_momentum_balance, storage_id * 1000 + i, JuMP.@constraint(gm.model, exp(beta) * rho_bottom^2 - rho_top^2 == (-resistance * phi_avg * abs(phi_avg)) * (exp(beta) - 1) / beta))
end
end

Expand Down
2 changes: 1 addition & 1 deletion src/core/transient_expression.jl
Original file line number Diff line number Diff line change
Expand Up @@ -159,7 +159,7 @@ function expression_compressor_power(gm::AbstractGasModel, nw::Int; report::Bool
m = (gm.ref[:it][gm_it_sym][:specific_heat_capacity_ratio] - 1) /
gm.ref[:it][gm_it_sym][:specific_heat_capacity_ratio]
W = 286.76 * gm.ref[:it][gm_it_sym][:temperature] / gm.ref[:it][gm_it_sym][:gas_specific_gravity] / m
var(gm, nw, :compressor_power_expr)[i] = JuMP.@NLexpression(gm.model, W * abs(f) * (alpha^m - 1.0))
var(gm, nw, :compressor_power_expr)[i] = JuMP.@expression(gm.model, W * abs(f) * (alpha^m - 1.0))
end

report && sol_component_value(
Expand Down
8 changes: 4 additions & 4 deletions src/form/dwp.jl
Original file line number Diff line number Diff line change
Expand Up @@ -236,8 +236,8 @@ function constraint_compressor_energy(gm::AbstractDWPModel, n::Int, k::Int, powe

M = abs(flow_max) * abs(ratio_max^2^m - 1) # working in the space of r^2

_add_constraint!(gm, n, :compressor_energy, k, JuMP.@NLconstraint(gm.model, f * (r^m - 1) <= power_max / work + ((1-y) * M)))
_add_constraint!(gm, n, :compressor_energy, k, JuMP.@NLconstraint(gm.model, -f * (r^m - 1) <= power_max / work + (y * M) ))
_add_constraint!(gm, n, :compressor_energy, k, JuMP.@constraint(gm.model, f * (r^m - 1) <= power_max / work + ((1-y) * M)))
_add_constraint!(gm, n, :compressor_energy, k, JuMP.@constraint(gm.model, -f * (r^m - 1) <= power_max / work + (y * M) ))
end


Expand All @@ -250,6 +250,6 @@ function constraint_compressor_energy_ne(gm::AbstractDWPModel, n::Int, k, power_
M = abs(flow_max) * abs(ratio_max^2^m - 1) # working in the space of r^2

# f is zero when the compressor is not built, so constraint is always true then
_add_constraint!(gm, n, :ne_compressor_energy, k, JuMP.@NLconstraint(gm.model, f * (r^m - 1) <= power_max / work + ((1-y) * M)))
_add_constraint!(gm, n, :ne_compressor_energy, k, JuMP.@NLconstraint(gm.model, -f * (r^m - 1) <= power_max / work + (y * M) ))
_add_constraint!(gm, n, :ne_compressor_energy, k, JuMP.@constraint(gm.model, f * (r^m - 1) <= power_max / work + ((1-y) * M)))
_add_constraint!(gm, n, :ne_compressor_energy, k, JuMP.@constraint(gm.model, -f * (r^m - 1) <= power_max / work + (y * M) ))
end
22 changes: 11 additions & 11 deletions src/form/wp.jl
Original file line number Diff line number Diff line change
Expand Up @@ -26,9 +26,9 @@ function constraint_pipe_weymouth(gm::AbstractWPModel, n::Int, k, i, j, f_min, f
f = var(gm, n, :f_pipe, k)

if w == 0.0
_add_constraint!(gm, n, :weymouth1, k, JuMP.@NLconstraint(gm.model, (pi - pj) == 0.0))
_add_constraint!(gm, n, :weymouth1, k, JuMP.@constraint(gm.model, (pi - pj) == 0.0))
else
_add_constraint!(gm, n, :weymouth1, k, JuMP.@NLconstraint(gm.model, (pi - pj) == (f * abs(f)) / w))
_add_constraint!(gm, n, :weymouth1, k, JuMP.@constraint(gm.model, (pi - pj) == (f * abs(f)) / w))
end
end

Expand All @@ -41,7 +41,7 @@ function constraint_inclined_pipe_pressure_drop(gm::AbstractWPModel, n::Int, k,
pj = var(gm, n, :psqr, j)
f = var(gm, n, :f_pipe, k)

_add_constraint!(gm, n, :inclined_pipe_pressure_drop, k, JuMP.@NLconstraint(gm.model, pj - exp(r_2)*pii== (exp(r_2) - 1)*r_1*f*abs(f)))
_add_constraint!(gm, n, :inclined_pipe_pressure_drop, k, JuMP.@constraint(gm.model, pj - exp(r_2)*pii== (exp(r_2) - 1)*r_1*f*abs(f)))
end


Expand All @@ -50,7 +50,7 @@ function constraint_resistor_darcy_weisbach(gm::AbstractWPModel, n::Int, k, i, j
f = var(gm, n, :f_resistor, k)
p_i, p_j = var(gm, n, :p, i), var(gm, n, :p, j)

_add_constraint!(gm, n, :darcy_weisbach_1, k, JuMP.@NLconstraint(gm.model, (1.0/w) * (p_i - p_j) == f * abs(f)))
_add_constraint!(gm, n, :darcy_weisbach_1, k, JuMP.@constraint(gm.model, (1.0/w) * (p_i - p_j) == f * abs(f)))
end


Expand Down Expand Up @@ -91,14 +91,14 @@ function constraint_pipe_weymouth_ne(gm::AbstractWPModel, n::Int, k, i, j, w, f_
# when z = 1, constraint is active
# when z = 0 f is also 0. Therefore, the big M we need is just the smallest and largest pressure difference that is possible
aux = JuMP.@variable(gm.model)
JuMP.@NLconstraint(gm.model, aux == f * abs(f))
JuMP.@constraint(gm.model, aux == f * abs(f))

if w == 0.0
_add_constraint!(gm, n, :weymouth_ne1, k, JuMP.@constraint(gm.model, pi - pj <= (1 - zp) * pd_max))
_add_constraint!(gm, n, :weymouth_ne2, k, JuMP.@constraint(gm.model, pi - pj >= (1 - zp) * pd_min))
else
_add_constraint!(gm, n, :weymouth_ne1, k, JuMP.@NLconstraint(gm.model, (pi - pj) <= aux / w + (1 - zp) * pd_max))
_add_constraint!(gm, n, :weymouth_ne2, k, JuMP.@NLconstraint(gm.model, (pi - pj) >= aux / w + (1 - zp) * pd_min))
_add_constraint!(gm, n, :weymouth_ne1, k, JuMP.@constraint(gm.model, (pi - pj) <= aux / w + (1 - zp) * pd_max))
_add_constraint!(gm, n, :weymouth_ne2, k, JuMP.@constraint(gm.model, (pi - pj) >= aux / w + (1 - zp) * pd_min))
end
end

Expand Down Expand Up @@ -243,7 +243,7 @@ function constraint_compressor_ratio_value(gm::AbstractWPModel, n::Int, k, i, j,
if type == 0
if (min_ratio <= 1.0 && max_ratio >= 1)
pk = get_compressor_pressure_aux(gm, n, k)
_add_constraint!(gm, n, :compressor_ratio_value1, k, JuMP.@NLconstraint(gm.model, pi*pj*r == pk*pj + pk*pi - pj*pi))
_add_constraint!(gm, n, :compressor_ratio_value1, k, JuMP.@constraint(gm.model, pi*pj*r == pk*pj + pk*pi - pj*pi))
else
y = get_compressor_y(gm, n, k)
_add_constraint!(gm, n, :compressor_ratio_value1, k, JuMP.@constraint(gm.model, r * pi <= pj + (1 - y) * i_pmax^2 * max_ratio^2))
Expand All @@ -268,7 +268,7 @@ function constraint_compressor_ratio_value_ne(gm::AbstractWPModel, n::Int, k, i,
if type == 0
if (min_ratio <= 1.0 && max_ratio >= 1)
pk = get_ne_compressor_pressure_aux(gm, n, k)
_add_constraint!(gm, n, :compressor_ratio_value1, k, JuMP.@NLconstraint(gm.model, pi*pj*r == pk*pj + pk*pi - pj*pi))
_add_constraint!(gm, n, :compressor_ratio_value1, k, JuMP.@constraint(gm.model, pi*pj*r == pk*pj + pk*pi - pj*pi))
else
y = get_ne_compressor_y_wp(gm, n, k)
_add_constraint!(gm, n, :compressor_ratio_value_ne1, k, JuMP.@constraint(gm.model, r * pi <= pj + (2 - y - z) * i_pmax^2 * max_ratio^2))
Expand All @@ -288,7 +288,7 @@ end
function constraint_compressor_energy(gm::AbstractWPModel, n::Int, k, power_max, m, work, flow_max, ratio_max)
r = var(gm, n, :rsqr, k)
f = var(gm, n, :f_compressor, k)
_add_constraint!(gm, n, :compressor_energy, k, JuMP.@NLconstraint(gm.model, abs(f) * (r^m - 1) <= power_max / work))
_add_constraint!(gm, n, :compressor_energy, k, JuMP.@constraint(gm.model, abs(f) * (r^m - 1) <= power_max / work))
end


Expand All @@ -297,5 +297,5 @@ function constraint_compressor_energy_ne(gm::AbstractWPModel, n::Int, k, power_m
r = var(gm, n, :rsqr_ne, k)
f = var(gm, n, :f_ne_compressor, k)
# when the compressor is not built, f = 0, so this is always true
_add_constraint!(gm, n, :ne_compressor_energy, k, JuMP.@NLconstraint(gm.model, abs(f) * (r^m - 1) <= power_max / work))
_add_constraint!(gm, n, :ne_compressor_energy, k, JuMP.@constraint(gm.model, abs(f) * (r^m - 1) <= power_max / work))
end

0 comments on commit 7a0dc25

Please sign in to comment.