Skip to content

Commit a10bcdc

Browse files
authored
Merge pull request lanl-ansi#263 from lanl-ansi/complementarity-constraints
Introduction Complementarity constraints for Weymouth Equations
2 parents 107adce + 5817ee3 commit a10bcdc

File tree

8 files changed

+290
-23
lines changed

8 files changed

+290
-23
lines changed

src/GasModels.jl

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -77,6 +77,8 @@ module GasModels
7777
include("form/wp.jl")
7878
include("form/dwp.jl")
7979
include("form/crdwp.jl")
80+
include("form/cwp.jl")
81+
8082

8183
include("prob/gf.jl")
8284
include("prob/ne.jl")

src/core/types.jl

Lines changed: 10 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -2,42 +2,44 @@
22
"DWP models"
33
abstract type AbstractDWPModel <: AbstractGasModel end
44

5-
65
"LRDWP models"
76
abstract type AbstractLRDWPModel <: AbstractGasModel end
87

9-
108
"CRDWP models"
119
abstract type AbstractCRDWPModel <: AbstractGasModel end
1210

13-
1411
"WP models"
1512
abstract type AbstractWPModel <: AbstractGasModel end
1613

14+
"CWP models"
15+
abstract type AbstractCWPModel <: AbstractWPModel end
1716

1817
"LRWP models"
1918
abstract type AbstractLRWPModel <: AbstractGasModel end
2019

21-
2220
"LRWP Model Type"
2321
mutable struct LRWPGasModel <: AbstractLRWPModel @gm_fields end
2422

25-
2623
"LRDWP Model Type"
2724
mutable struct LRDWPGasModel <: AbstractLRDWPModel @gm_fields end
2825

29-
3026
"WP Model Type"
3127
mutable struct WPGasModel <: AbstractWPModel @gm_fields end
3228

33-
3429
"CRDWP Model Type"
3530
mutable struct CRDWPGasModel <: AbstractCRDWPModel @gm_fields end
3631

37-
3832
"DWP Model Type"
3933
mutable struct DWPGasModel <: AbstractDWPModel @gm_fields end
4034

35+
"""
36+
ComplementarityWPModel <: AbstractWPModel
37+
38+
Custom GasModels.jl model type for a Weymouth formulation that uses two
39+
nonnegative flow variables (f⁺ and f⁻) per pipe, plus a complementarity
40+
constraint f⁺·f⁻ = 0.
41+
"""
42+
mutable struct CWPGasModel <: AbstractCWPModel @gm_fields end
4143

4244
"Union of MI Models"
4345
AbstractMIModels = Union{AbstractCRDWPModel,AbstractDWPModel,AbstractLRDWPModel}

src/form/cwp.jl

Lines changed: 190 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,190 @@
1+
2+
"Variables needed for modeling flow in CWP models"
3+
function variable_flow(gm::AbstractCWPModel, nw::Int = nw_id_default; bounded::Bool = true, report::Bool = true)
4+
variable_mass_flow(gm, nw; bounded = bounded, report = report)
5+
variable_pipe_mass_flow_slack(gm, nw; bounded, report)
6+
end
7+
8+
9+
"Variables needed for modeling flow in CWP models"
10+
function variable_flow_ne(gm::AbstractCWPModel, nw::Int = nw_id_default; bounded::Bool = true, report::Bool = true)
11+
variable_mass_flow_ne(gm, nw; bounded = bounded, report = report)
12+
variable_pipe_mass_flow_slack_ne(gm, nw; bounded, report)
13+
end
14+
15+
16+
"""
17+
variable_pipe_mass_flow(gm::CWPGasModel, data, cdata)
18+
19+
Overrides the default pipe flow variable creation to define complementarity-based
20+
flow variables f⁺[p] and f⁻[p], along with a net flow variable f[p] = f⁺[p] - f⁻[p].
21+
"""
22+
function variable_pipe_mass_flow_slack(gm::AbstractCWPModel, nw::Int=nw_id_default; bounded::Bool=true, report::Bool=true)
23+
# 1) Define f⁺ variables (nonnegative)
24+
f_plus_pipe = var(gm, nw)[:f_plus_pipe] = JuMP.@variable(
25+
gm.model,
26+
[i in ids(gm, nw, :pipe)],
27+
lower_bound = 0,
28+
base_name = "$(nw)_f_plus",
29+
start = comp_start_value(ref(gm, nw, :pipe), i, "f_plus_start", 0)
30+
)
31+
32+
# 2) Define f⁻ variables (nonnegative)
33+
f_minus_pipe = var(gm, nw)[:f_minus_pipe] = JuMP.@variable(
34+
gm.model,
35+
[i in ids(gm, nw, :pipe)],
36+
lower_bound = 0,
37+
base_name = "$(nw)_f_minus",
38+
start = comp_start_value(ref(gm, nw, :pipe), i, "f_minus_start", 0)
39+
)
40+
41+
# 2) Optionally set bounds from pipe data
42+
if bounded
43+
for (i, pipe_dict) in ref(gm, nw, :pipe)
44+
# handle flow_min
45+
if haskey(pipe_dict, "flow_min")
46+
flow_min = pipe_dict["flow_min"]
47+
if flow_min > 0
48+
# pipe demands a strictly positive lower bound on "forward" flow
49+
JuMP.set_lower_bound(f_plus_pipe[i], flow_min)
50+
elseif flow_min <= 0
51+
# pipe demands a strictly negative flow_min => set a bound on "reverse" flow
52+
JuMP.set_upper_bound(f_minus_pipe[i], -flow_min)
53+
end
54+
end
55+
56+
# handle flow_max
57+
if haskey(pipe_dict, "flow_max")
58+
flow_max = pipe_dict["flow_max"]
59+
if flow_max < 0
60+
# maximum is negative => limit the reverse flow variable
61+
JuMP.set_lower_bound(f_minus_pipe[i], -flow_max)
62+
elseif flow_max >= 0
63+
# standard positive flow bound
64+
JuMP.set_upper_bound(f_plus_pipe[i], flow_max)
65+
end
66+
end
67+
end
68+
end
69+
70+
# 4) Optionally register these variables for solution reporting
71+
if report
72+
sol_component_value(gm, nw, :pipe, :f_plus_pipe, ids(gm, nw, :pipe), f_plus_pipe)
73+
sol_component_value(gm, nw, :pipe, :f_minus_pipe, ids(gm, nw, :pipe), f_minus_pipe)
74+
end
75+
end
76+
77+
78+
"""
79+
variable_pipe_mass_flow(gm::CWPGasModel, data, cdata)
80+
81+
Overrides the default pipe flow variable creation to define complementarity-based
82+
flow variables f⁺[p] and f⁻[p], along with a net flow variable f[p] = f⁺[p] - f⁻[p].
83+
"""
84+
function variable_pipe_mass_flow_slack_ne(gm::AbstractCWPModel, nw::Int=nw_id_default; bounded::Bool=true, report::Bool=true)
85+
# 1) Define f⁺ variables (nonnegative)
86+
f_plus_pipe = var(gm, nw)[:f_plus_ne_pipe] = JuMP.@variable(
87+
gm.model,
88+
[i in ids(gm, nw, :ne_pipe)],
89+
lower_bound = 0,
90+
base_name = "$(nw)_f_plus",
91+
start = comp_start_value(ref(gm, nw, :ne_pipe), i, "f_plus_start", 0)
92+
)
93+
94+
# 2) Define f⁻ variables (nonnegative)
95+
f_minus_pipe = var(gm, nw)[:f_minus_ne_pipe] = JuMP.@variable(
96+
gm.model,
97+
[i in ids(gm, nw, :ne_pipe)],
98+
lower_bound = 0,
99+
base_name = "$(nw)_f_minus",
100+
start = comp_start_value(ref(gm, nw, :ne_pipe), i, "f_minus_start", 0)
101+
)
102+
103+
# 2) Optionally set bounds from pipe data
104+
if bounded
105+
for (i, pipe_dict) in ref(gm, nw, :ne_pipe)
106+
# handle flow_min
107+
if haskey(pipe_dict, "flow_min")
108+
flow_min = pipe_dict["flow_min"]
109+
if flow_min <= 0
110+
# pipe demands a strictly negative flow_min => set a bound on "reverse" flow
111+
JuMP.set_upper_bound(f_minus_pipe[i], -flow_min)
112+
end
113+
end
114+
115+
# handle flow_max
116+
if haskey(pipe_dict, "flow_max")
117+
flow_max = pipe_dict["flow_max"]
118+
if flow_max >= 0
119+
# standard positive flow bound
120+
JuMP.set_upper_bound(f_plus_pipe[i], flow_max)
121+
end
122+
end
123+
end
124+
end
125+
126+
# 4) Optionally register these variables for solution reporting
127+
if report
128+
sol_component_value(gm, nw, :ne_pipe, :f_plus_ne_pipe, ids(gm, nw, :ne_pipe), f_plus_pipe)
129+
sol_component_value(gm, nw, :ne_pipe, :f_minus_ne_pipe, ids(gm, nw, :ne_pipe), f_minus_pipe)
130+
end
131+
end
132+
133+
"Weymouth equation with absolute value"
134+
function constraint_pipe_weymouth(gm::AbstractCWPModel, n::Int, k, i, j, f_min, f_max, w, pd_min, pd_max)
135+
# Retrieve squared pressures at the 'from' and 'to' nodes
136+
pi = var(gm, n, :psqr, i)
137+
pj = var(gm, n, :psqr, j)
138+
139+
# Retrieve f for this pipe
140+
f_plus = var(gm, n, :f_plus_pipe, k)
141+
f_minus = var(gm, n, :f_minus_pipe, k)
142+
f = var(gm, n, :f_pipe, k)
143+
144+
if w == 0.0
145+
# Degenerate case
146+
_add_constraint!(gm, n, :weymouth1, k, JuMP.@constraint(gm.model, pi - pj == 0.0))
147+
else
148+
_add_constraint!(gm, n, :weymouth1, k, JuMP.@constraint(gm.model, pi - pj == (f_plus^2 - f_minus^2) / w))
149+
end
150+
151+
# Complementarity (no simultaneous forward and reverse flow):
152+
_add_constraint!(gm, n, :complementarity_weymouth, k, JuMP.@constraint(gm.model, f_plus * f_minus == 0))
153+
154+
#Relate the net flow variable f to f_plus and f_minus:
155+
_add_constraint!(gm, n, :weymouth_flow_relation, k, JuMP.@constraint(gm.model, f == f_plus - f_minus))
156+
end
157+
158+
"Weymouth equation for an expansion pipe"
159+
function constraint_pipe_weymouth_ne(gm::AbstractCWPModel, n::Int, k, i, j, w, f_min, f_max, pd_min, pd_max)
160+
# Retrieve squared pressures at the 'from' and 'to' nodes
161+
pi = var(gm, n, :psqr, i)
162+
pj = var(gm, n, :psqr, j)
163+
164+
# Retrieve f for this pipe
165+
f_plus = var(gm, n, :f_plus_ne_pipe, k)
166+
f_minus = var(gm, n, :f_minus_ne_pipe, k)
167+
f = var(gm, n, :f_ne_pipe, k)
168+
zp = var(gm, n, :zp, k)
169+
170+
# when z = 1, constraint is active
171+
# when z = 0 f is also 0. Therefore, the big M we need is just the smallest and largest pressure difference that is possible
172+
aux = JuMP.@variable(gm.model)
173+
JuMP.@constraint(gm.model, aux == f_plus^2 - f_minus^2)
174+
175+
if w == 0.0
176+
# Degenerate case
177+
_add_constraint!(gm, n, :weymouth_ne1, k, JuMP.@constraint(gm.model, pi - pj <= (1 - zp) * pd_max))
178+
_add_constraint!(gm, n, :weymouth_ne2, k, JuMP.@constraint(gm.model, pi - pj >= (1 - zp) * pd_min))
179+
else
180+
_add_constraint!(gm, n, :weymouth_ne1, k, JuMP.@constraint(gm.model, (pi - pj) <= aux / w + (1 - zp) * pd_max))
181+
_add_constraint!(gm, n, :weymouth_ne2, k, JuMP.@constraint(gm.model, (pi - pj) >= aux / w + (1 - zp) * pd_min))
182+
end
183+
184+
# Complementarity (no simultaneous forward and reverse flow):
185+
_add_constraint!(gm, n, :complementarity_weymouth, k, JuMP.@constraint(gm.model, f_plus * f_minus == 0))
186+
187+
#Relate the net flow variable f to f_plus and f_minus:
188+
_add_constraint!(gm, n, :weymouth_flow_relation, k, JuMP.@constraint(gm.model, f == f_plus - f_minus))
189+
end
190+

test/gf.jl

Lines changed: 10 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -37,6 +37,16 @@
3737
@test result["termination_status"] in [LOCALLY_SOLVED, ALMOST_LOCALLY_SOLVED, OPTIMAL, :Suboptimal]
3838
end
3939

40+
@testset "test cwp gf" begin
41+
@info "Testing cwp gf"
42+
result = run_gf("../test/data/matgas/case-6-gf.m", CWPGasModel, nlp_solver)
43+
@test result["termination_status"] in [LOCALLY_SOLVED, ALMOST_LOCALLY_SOLVED, OPTIMAL, :Suboptimal]
44+
@test isapprox(result["objective"], 0; atol = 1e-6)
45+
46+
result = run_gf("../test/data/gaslib/GasLib-Integration.zip", WPGasModel, minlp_solver)
47+
@test result["termination_status"] in [LOCALLY_SOLVED, ALMOST_LOCALLY_SOLVED, OPTIMAL, :Suboptimal]
48+
end
49+
4050
@testset "test dwp gf" begin
4151
@info "Testing dwp gf"
4252
result = run_gf("../test/data/matgas/case-6-gf.m", DWPGasModel, minlp_solver)

test/ls.jl

Lines changed: 31 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -1,65 +1,81 @@
11
@testset "test ls" begin
22
#Check the second order cone model on load shedding
33
@testset "test crdwp ls" begin
4-
@info "Testing gaslib crdwp ls"
4+
@info "Testing case 6 crdwp ls"
55
result = run_ls("../test/data/matgas/case-6-ls.m", CRDWPGasModel, misocp_solver)
66
@test result["termination_status"] in [LOCALLY_SOLVED, ALMOST_LOCALLY_SOLVED, OPTIMAL, :Suboptimal]
77
@test isapprox(result["objective"] * result["solution"]["base_flow"], 5.0; atol = 1e-1)
88
end
99

1010
#Check the second order cone model on load shedding with priorities
11-
@testset "test gaslib 40 crdwp ls priority" begin
12-
@info "Testing gaslib crdwp ls priority gaslib 40"
11+
@testset "test case 6 crdwp ls priority" begin
12+
@info "Testing case 6 crdwp ls priority gaslib 40"
1313
result = run_ls("../test/data/matgas/case-6-ls-priority.m", CRDWPGasModel, misocp_solver)
1414
@test result["termination_status"] in [LOCALLY_SOLVED, ALMOST_LOCALLY_SOLVED, OPTIMAL, :Suboptimal]
1515
@test isapprox(result["objective"] * result["solution"]["base_flow"], 4.5; atol = 1e-1)
1616
end
1717

1818
#Check the lrdwp model on load shedding
19-
@testset "test gaslib 40 lrdwp ls" begin
20-
@info "Testing gaslib lrdwp ls gaslib 40"
19+
@testset "test case 6 lrdwp ls" begin
20+
@info "Testing case 6 lrdwp ls gaslib 40"
2121
result = run_ls("../test/data/matgas/case-6-ls.m", LRDWPGasModel, mip_solver)
2222
@test result["termination_status"] in [LOCALLY_SOLVED, ALMOST_LOCALLY_SOLVED, OPTIMAL, :Suboptimal]
2323
@test isapprox(result["objective"] * result["solution"]["base_flow"], 5.0; atol = 1e-1)
2424
end
2525

2626
#Check the lrdwp model on load shedding with priorities
27-
@testset "test gaslib 40 lrdwp ls priority" begin
28-
@info "Testing gaslib lrdwp ls priority gaslib 40"
27+
@testset "test case 6 lrdwp ls priority" begin
28+
@info "Testing case 6 lrdwp ls priority gaslib 40"
2929
result = run_ls("../test/data/matgas/case-6-ls-priority.m", LRDWPGasModel, mip_solver)
3030
@test result["termination_status"] in [LOCALLY_SOLVED, ALMOST_LOCALLY_SOLVED, OPTIMAL, :Suboptimal]
3131
@test isapprox(result["objective"] * result["solution"]["base_flow"], 4.5; atol = 1e-1)
3232
end
3333

3434
#Check the lrwp model on load shedding
35-
@testset "test gaslib 40 lrwp ls" begin
36-
@info "Testing gaslib lrwp ls gaslib 40"
35+
@testset "test case 6 lrwp ls" begin
36+
@info "Testing case 6 lrwp ls gaslib 40"
3737
result = run_ls("../test/data/matgas/case-6-ls.m", LRWPGasModel, lp_solver)
3838
@test result["termination_status"] in [LOCALLY_SOLVED, ALMOST_LOCALLY_SOLVED, OPTIMAL, :Suboptimal]
3939
@test isapprox(result["objective"] * result["solution"]["base_flow"], 5.0; atol = 1e-1)
4040
end
4141

4242
#Check the lrwp model on load shedding with priorities
43-
@testset "test gaslib 40 lrwp ls priority" begin
44-
@info "Testing gaslib lrwp ls priority gaslib 40"
43+
@testset "test case 6 lrwp ls priority" begin
44+
@info "Testing case 6 lrwp ls priority gaslib 40"
4545
result = run_ls("../test/data/matgas/case-6-ls-priority.m", LRWPGasModel, lp_solver)
4646
@test result["termination_status"] in [LOCALLY_SOLVED, ALMOST_LOCALLY_SOLVED, OPTIMAL, :Suboptimal]
4747
@test isapprox(result["objective"] * result["solution"]["base_flow"], 4.5; atol = 1e-1)
4848
end
4949

5050
#Check the wp model on load shedding
51-
@testset "test gaslib 40 wp ls" begin
52-
@info "Testing gaslib wp ls gaslib 40"
51+
@testset "test case 6 wp ls" begin
52+
@info "Testing case 6 wp ls case"
5353
result = run_ls("../test/data/matgas/case-6-ls.m", WPGasModel, nlp_solver)
5454
@test result["termination_status"] in [LOCALLY_SOLVED, ALMOST_LOCALLY_SOLVED, OPTIMAL, :Suboptimal]
5555
@test isapprox(result["objective"] * result["solution"]["base_flow"], 5.0; atol = 1e-1)
5656
end
5757

5858
#Check the wp model on load shedding with priorities
59-
@testset "test gaslib 40 wp priority" begin
60-
@info "Testing gaslib wp ls priority gaslib 40"
59+
@testset "test case 6 wp priority" begin
60+
@info "Testing case 6 ls priority case 6"
6161
result = run_ls("../test/data/matgas/case-6-ls-priority.m", WPGasModel, nlp_solver)
6262
@test result["termination_status"] in [LOCALLY_SOLVED, ALMOST_LOCALLY_SOLVED, OPTIMAL, :Suboptimal]
6363
@test isapprox(result["objective"] * result["solution"]["base_flow"], 4.5; atol = 1e-1)
6464
end
65+
66+
#Check the cwp model on load shedding
67+
@testset "test case 6 cwp ls" begin
68+
@info "Testing case 6 cwp ls gaslib 40"
69+
result = run_ls("../test/data/matgas/case-6-ls.m", CWPGasModel, nlp_solver)
70+
@test result["termination_status"] in [LOCALLY_SOLVED, ALMOST_LOCALLY_SOLVED, OPTIMAL, :Suboptimal]
71+
@test isapprox(result["objective"] * result["solution"]["base_flow"], 5.0; atol = 1e-1)
72+
end
73+
74+
#Check the cwp model on load shedding with priorities
75+
@testset "test case 6 40 cwp priority" begin
76+
@info "Testing case 6 cwp ls priority case 6"
77+
result = run_ls("../test/data/matgas/case-6-ls-priority.m", CWPGasModel, nlp_solver)
78+
@test result["termination_status"] in [LOCALLY_SOLVED, ALMOST_LOCALLY_SOLVED, OPTIMAL, :Suboptimal]
79+
@test isapprox(result["objective"] * result["solution"]["base_flow"], 4.5; atol = 1e-1)
80+
end
6581
end

test/ne.jl

Lines changed: 10 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -36,4 +36,14 @@
3636
@test result["termination_status"] in [ALMOST_LOCALLY_SOLVED, OPTIMAL, :Suboptimal]
3737
end
3838
end
39+
40+
@testset "test cwp ne" begin
41+
@info "Testing cwp ne"
42+
result = run_ne("../test/data/matgas/case-6-ne.m", CWPGasModel, minlp_solver)
43+
if result["termination_status"] == LOCALLY_SOLVED
44+
@test isapprox(result["objective"], 1476; atol = 1e-1)
45+
else # CI compat for windows on Julia v1.6, 01/29/24
46+
@test result["termination_status"] in [ALMOST_LOCALLY_SOLVED, OPTIMAL, :Suboptimal]
47+
end
48+
end
3949
end

test/nels.jl

Lines changed: 7 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -33,4 +33,11 @@
3333
@test result["termination_status"] in [LOCALLY_SOLVED, ALMOST_LOCALLY_SOLVED, OPTIMAL, :Suboptimal]
3434
@test isapprox(result["objective"] * result["solution"]["base_flow"], 1031.60; atol = 1e-1)
3535
end
36+
37+
@testset "test cwp nels" begin
38+
@info "Testing cwp nels"
39+
result = run_nels("../test/data/matgas/case-6-nels.m", CWPGasModel, minlp_solver)
40+
@test result["termination_status"] in [LOCALLY_SOLVED, ALMOST_LOCALLY_SOLVED, OPTIMAL, :Suboptimal]
41+
@test isapprox(result["objective"] * result["solution"]["base_flow"], 1031.60; atol = 1e-1)
42+
end
3643
end

0 commit comments

Comments
 (0)