From c268bfe637292a5c0337aa4c85e10e76948f001f Mon Sep 17 00:00:00 2001 From: Sai Krishna Kanth Hari Date: Fri, 9 Sep 2022 12:18:47 -0700 Subject: [PATCH 1/7] Correcting reservoir pressure update Pressure is supposed to decrease with withdrawal of gas --- src/core/constraint_transient.jl | 2 +- test/transient.jl | 16 ++++++++-------- 2 files changed, 9 insertions(+), 9 deletions(-) diff --git a/src/core/constraint_transient.jl b/src/core/constraint_transient.jl index 9a582a4d..b55822da 100644 --- a/src/core/constraint_transient.jl +++ b/src/core/constraint_transient.jl @@ -131,7 +131,7 @@ function constraint_storage_reservoir_physics_simplified( f = 0.5 * (var(gm, nw, :storage_flow, storage_id) + var(gm, nw + 1, :storage_flow, storage_id)) end - GasModels._add_constraint!(gm, nw, :reservoir_physics, storage_id, JuMP.@constraint(gm.model, volume * rho_dot == f)) + GasModels._add_constraint!(gm, nw, :reservoir_physics, storage_id, JuMP.@constraint(gm.model, volume * rho_dot == -f)) end "Constraint: flow bounds imposed by injection/withdrawal well" diff --git a/test/transient.jl b/test/transient.jl index 0570ba28..be9cac68 100644 --- a/test/transient.jl +++ b/test/transient.jl @@ -41,25 +41,25 @@ end end @testset "transient (steady state replicate) case with storage" begin - mn_data = parse_files("../test/data/matgas/case-6-storage.m", "../test/data/transient/time-series-case-6b.csv", spatial_discretization = 1e4, additional_time = 3600.0) + mn_data = parse_files("../test/data/matgas/case-6-storage.m", "../test/data/transient/time-series-case-6b.csv", spatial_discretization = 1e4, additional_time = 7200.0) result = run_transient_ogf(mn_data, WPGasModel, nlp_solver) @test result["termination_status"] == LOCALLY_SOLVED make_si_units!(result["solution"]) - @test isapprox(result["solution"]["nw"]["1"]["receipt"]["1"]["injection"], 41, atol = 1) - @test isapprox(result["solution"]["nw"]["2"]["receipt"]["1"]["injection"], 41, atol = 1) - @test isapprox(result["solution"]["nw"]["1"]["storage"]["1"]["storage_flow"], 88, atol = 1) - @test isapprox(result["solution"]["nw"]["2"]["storage"]["1"]["storage_flow"], 88, atol = 1) + @test isapprox(result["solution"]["nw"]["1"]["receipt"]["1"]["injection"], 43.41, atol = 1) + @test isapprox(result["solution"]["nw"]["2"]["receipt"]["1"]["injection"], 43.02, atol = 1) + @test isapprox(result["solution"]["nw"]["1"]["storage"]["1"]["storage_flow"], 86.45, atol = 1) + @test isapprox(result["solution"]["nw"]["2"]["storage"]["1"]["storage_flow"], 87.96, atol = 1) end @testset "transient time-periodic withdrawal case with storage" begin mn_data = parse_files("../test/data/matgas/case-6-storage.m", "../test/data/transient/time-series-case-6a.csv", spatial_discretization = 1e4, additional_time = 3600.0) result = run_transient_ogf(mn_data, WPGasModel, nlp_solver) - @test isapprox(result["objective"], -9282.09, atol = 1) + @test isapprox(result["objective"], -9220.2, atol = 1) make_si_units!(result["solution"]) @test isapprox(result["solution"]["nw"]["1"]["receipt"]["1"]["injection"], 0.0, atol = 1e-1) @test isapprox(result["solution"]["nw"]["2"]["receipt"]["1"]["injection"], 0.0, atol = 1e-1) - @test isapprox(result["solution"]["nw"]["1"]["storage"]["1"]["storage_flow"], 88, atol = 1) - @test isapprox(result["solution"]["nw"]["2"]["storage"]["1"]["storage_flow"], 88, atol = 1) + @test isapprox(result["solution"]["nw"]["1"]["storage"]["1"]["storage_flow"], 86.51, atol = 1) + @test isapprox(result["solution"]["nw"]["2"]["storage"]["1"]["storage_flow"], 87.96, atol = 1) end @testset "transient (steady state replicate) case with elevation" begin From 837d7fe3674ad23351ad2d2a586236176e5a0bf0 Mon Sep 17 00:00:00 2001 From: Sai Krishna Kanth Hari Date: Mon, 27 Nov 2023 21:11:31 -0600 Subject: [PATCH 2/7] lrwp for weymouth and new objective for compressor power Added lrwp version of the pipe weymouth constraint. Created a new ogf problem by replacing the compressor power objective with the difference in the squared pressures at the compressor end points --- Project.toml | 4 ++ src/GasModels.jl | 2 + src/core/objective.jl | 65 +++++++++++++++++++++++++ src/form/crdwp.jl | 3 ++ src/form/dwp.jl | 3 ++ src/form/lrdwp.jl | 4 +- src/form/lrwp.jl | 29 +++++++++++- src/form/wp.jl | 4 +- src/prob/gf.jl | 1 + src/prob/ls.jl | 1 + src/prob/ne.jl | 1 + src/prob/nels.jl | 1 + src/prob/ogf.jl | 1 + src/prob/ogf_new.jl | 107 ++++++++++++++++++++++++++++++++++++++++++ test/ogf.jl | 11 +++++ 15 files changed, 233 insertions(+), 4 deletions(-) create mode 100644 src/prob/ogf_new.jl diff --git a/Project.toml b/Project.toml index c3b31174..6edc5508 100644 --- a/Project.toml +++ b/Project.toml @@ -15,6 +15,9 @@ 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" @@ -29,6 +32,7 @@ XMLDict = "~0.4" ZipFile = "~0.8, ~0.9" julia = "^1" + [extras] HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b" MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" diff --git a/src/GasModels.jl b/src/GasModels.jl index 484e82f8..621a9c4c 100644 --- a/src/GasModels.jl +++ b/src/GasModels.jl @@ -12,6 +12,8 @@ module GasModels using Dates using Dierckx + using PolyhedralRelaxations + # Create our module level logger (this will get precompiled) diff --git a/src/core/objective.jl b/src/core/objective.jl index 7e5026cd..bcabfe99 100644 --- a/src/core/objective.jl +++ b/src/core/objective.jl @@ -200,3 +200,68 @@ function objective_min_transient_compressor_power(gm::AbstractGasModel, time_poi return JuMP.@objective(gm.model, Min, compressor_power_expression) end + + + +"function for minimizing new economic costs: ``\\min \\sum_{j \\in {\\cal D}} \\kappa_j \\boldsymbol{d}_j - \\sum_{j \\in {\\cal T}} \\kappa_j \\boldsymbol{\\tau}_j - \\sum_{j \\in {\\cal R}} \\kappa_j \\boldsymbol{r}_j - + \\sum_{k \\in {\\cal C}} \\boldsymbol{p^2}_{jk} - \\boldsymbol{p^2}_{ik} ``" +function objective_min_new_economic_costs(gm::AbstractGasModel, nws = [nw_id_default]) + r = Dict(n => var(gm, n, :rsqr) for n in nws) + f = Dict(n => var(gm, n, :f_compressor) for n in nws) + p2 = Dict(n => var(gm, n, :psqr) for n in nws) + fl = Dict(n => var(gm, n, :fl) for n in nws) + fg = Dict(n => var(gm, n, :fg) for n in nws) + ft = Dict(n => var(gm, n, :ft) for n in nws) + gamma = get_specific_heat_capacity_ratio(gm.data) + m = ((gamma - 1) / gamma) / 2 + T = get_temperature(gm.data) + G = get_gas_specific_gravity(gm.data) + work = _calc_compressor_work(gamma, G, T) + base_flow = get_base_flow(gm.data) + + load_set = Dict( + n => keys(Dict( + x for x in ref(gm, n, :delivery) if x.second["is_dispatchable"] == 1 + )) for n in nws + ) + transfer_set = Dict( + n => keys(Dict( + x for x in ref(gm, n, :transfer) if x.second["is_dispatchable"] == 1 + )) for n in nws + ) + prod_set = Dict( + n => keys(Dict( + x for x in ref(gm, n, :receipt) if x.second["is_dispatchable"] == 1 + )) for n in nws + ) + load_prices = Dict( + n => Dict( + i => get(ref(gm, n, :delivery, i), "bid_price", 1.0) for i in load_set[n] + ) for n in nws + ) + prod_prices = Dict( + n => Dict( + i => get(ref(gm, n, :receipt, i), "offer_price", 1.0) for i in prod_set[n] + ) for n in nws + ) + transfer_prices = Dict( + n => Dict( + i => get(ref(gm, n, :transfer, i), "bid_price", 1.0) for i in transfer_set[n] + ) for n in nws + ) + + economic_weighting = get_economic_weighting(gm.data) + + # 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.@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]) + + economic_weighting * sum(prod_prices[n][i] * fg[n][i] for i in prod_set[n]) + + (1.0 - economic_weighting) * + sum( (p2[n][compressor["to_junction"]] - p2[n][compressor["fr_junction"]]) for (i, compressor) in ref(gm, n, :compressor)) + for n in nws + )) + return JuMP.@objective(gm.model, Min, z) +end diff --git a/src/form/crdwp.jl b/src/form/crdwp.jl index 43f7dc17..156a7afa 100644 --- a/src/form/crdwp.jl +++ b/src/form/crdwp.jl @@ -1,4 +1,7 @@ # Define CRDWP implementations of Gas Models +function variable_form_specific(gm::AbstractCRDWPModel, nw::Int=nw_id_default; bounded::Bool=true, report::Bool=true) + # NONE +end "Variables needed for modeling flow in MI models" function variable_flow(gm::AbstractCRDWPModel, nw::Int = nw_id_default; bounded::Bool = true, report::Bool = true) diff --git a/src/form/dwp.jl b/src/form/dwp.jl index 12b2caa9..aec77b12 100644 --- a/src/form/dwp.jl +++ b/src/form/dwp.jl @@ -1,4 +1,7 @@ # Define DWP implementations of Gas Models +function variable_form_specific(gm::AbstractDWPModel, nw::Int=nw_id_default; bounded::Bool=true, report::Bool=true) + # NONE +end "Variables needed for modeling flow in MI models" function variable_flow(gm::AbstractDWPModel, n::Int = nw_id_default; bounded::Bool = true, report::Bool = true) diff --git a/src/form/lrdwp.jl b/src/form/lrdwp.jl index 258dffed..7de4f7c2 100644 --- a/src/form/lrdwp.jl +++ b/src/form/lrdwp.jl @@ -1,5 +1,7 @@ # Define LRDWP implementations of Gas Models - +function variable_form_specific(gm::AbstractLRDWPModel, nw::Int=nw_id_default; bounded::Bool=true, report::Bool=true) + # NONE +end "Variables needed for modeling flow in LRDWP models" function variable_flow(gm::AbstractLRDWPModel, n::Int = nw_id_default; bounded::Bool = true, report::Bool = true) diff --git a/src/form/lrwp.jl b/src/form/lrwp.jl index 8433afd1..ea274928 100644 --- a/src/form/lrwp.jl +++ b/src/form/lrwp.jl @@ -1,12 +1,37 @@ # Define LRWP implementations of Gas Models - +function variable_form_specific(gm::AbstractLRWPModel, nw::Int=nw_id_default; bounded::Bool=true, report::Bool=true) + # pipe f|f\ relaxation + fmf_l_pipe = var(gm, nw)[:fmf_l_pipe] = JuMP.@variable(gm.model, + [i in ids(gm, nw, :pipe)], + base_name="$(nw)_fmf_l") + report && sol_component_value(gm, nw, :pipe, :fmf_l, ids(gm, nw, :pipe), fmf_l_pipe) +end ###################################################################################################### ## Constraints ###################################################################################################### "Constraint: Weymouth equation--not applicable for LRWP models" function constraint_pipe_weymouth(gm::AbstractLRWPModel, n::Int, k, i, j, f_min, f_max, w, pd_min, pd_max) - # TODO Linear convex hull of the weymouth equations in wp.jl + pi = var(gm, n, :psqr, i) + pj = var(gm, n, :psqr, j) + f = var(gm, n, :f_pipe, k) + fmf_l = var(gm, n, :fmf_l_pipe, k) + if w == 0.0 + _add_constraint!(gm, n, :weymouth1, k, JuMP.@constraint(gm.model, (pi - pj) == 0.0)) + elseif f_min == f_max + _add_constraint!(gm, n, :weymouth1, k, JuMP.@constraint(gm.model, w * (pi - pj) == f_min)) + else + _add_constraint!(gm, n, :weymouth1, k, JuMP.@constraint(gm.model, w * (pi - pj) == fmf_l)) + # fmf_f incorporates the univariate relaxation for f*(abs(f)) + if(f_min<0 && f_max>0) + partition = [f_min, 0, f_max] + # partition = [f_min,3*f_min/4,f_min/2,f_min/4,0,f_max/4,f_max/2,3*f_max/4,f_max] + else + partition = [f_min, f_max] + end + construct_univariate_relaxation!(gm.model, a -> a*(abs(a)), f, fmf_l, partition, false) + end + end diff --git a/src/form/wp.jl b/src/form/wp.jl index cbfe5768..6b134f9e 100644 --- a/src/form/wp.jl +++ b/src/form/wp.jl @@ -3,7 +3,9 @@ ############################################################################################################# ## Constraints for modeling flow across a pipe ############################################################################################################ - +function variable_form_specific(gm::AbstractWPModel, nw::Int=nw_id_default; bounded::Bool=true, report::Bool=true) + # NONE +end "Constraint: Constraints which define pressure drop across a resistor" function constraint_resistor_pressure(gm::AbstractWPModel, n::Int, k::Int, i::Int, j::Int, pd_min::Float64, pd_max::Float64) diff --git a/src/prob/gf.jl b/src/prob/gf.jl index f79c7377..45e0fae5 100644 --- a/src/prob/gf.jl +++ b/src/prob/gf.jl @@ -50,6 +50,7 @@ function build_gf(gm::AbstractGasModel) variable_transfer_mass_flow(gm) variable_compressor_ratio_sqr(gm; compressors = bounded_compressors) variable_storage(gm) + variable_form_specific(gm) for (i, junction) in ref(gm, :junction) constraint_mass_flow_balance(gm, i) diff --git a/src/prob/ls.jl b/src/prob/ls.jl index bdc67888..70a3290b 100644 --- a/src/prob/ls.jl +++ b/src/prob/ls.jl @@ -38,6 +38,7 @@ function build_ls(gm::AbstractGasModel) variable_transfer_mass_flow(gm) variable_compressor_ratio_sqr(gm; compressors = bounded_compressors) variable_storage(gm) + variable_form_specific(gm) objective_max_load(gm) diff --git a/src/prob/ne.jl b/src/prob/ne.jl index c6f100a3..5abd0177 100644 --- a/src/prob/ne.jl +++ b/src/prob/ne.jl @@ -54,6 +54,7 @@ function build_ne(gm::AbstractGasModel) variable_compressor_ratio_sqr(gm; compressors = bounded_compressors) variable_compressor_ratio_sqr_ne(gm; compressors = bounded_compressors_ne) variable_storage(gm) + variable_form_specific(gm) # expansion cost objective objective_min_ne_cost(gm) diff --git a/src/prob/nels.jl b/src/prob/nels.jl index 2f876e2f..00512ccc 100644 --- a/src/prob/nels.jl +++ b/src/prob/nels.jl @@ -51,6 +51,7 @@ function build_nels(gm::AbstractGasModel) variable_compressor_ratio_sqr(gm; compressors = bounded_compressors) variable_compressor_ratio_sqr_ne(gm; compressors = bounded_compressors_ne) variable_storage(gm) + variable_form_specific(gm) # expansion variables variable_pipe_ne(gm) diff --git a/src/prob/ogf.jl b/src/prob/ogf.jl index b4c93451..4270c416 100644 --- a/src/prob/ogf.jl +++ b/src/prob/ogf.jl @@ -50,6 +50,7 @@ function build_ogf(gm::AbstractGasModel) variable_transfer_mass_flow(gm) variable_compressor_ratio_sqr(gm) variable_storage(gm) + variable_form_specific(gm) objective_min_economic_costs(gm) diff --git a/src/prob/ogf_new.jl b/src/prob/ogf_new.jl new file mode 100644 index 00000000..c0c20d1d --- /dev/null +++ b/src/prob/ogf_new.jl @@ -0,0 +1,107 @@ +# Definitions for running the new optimal gas flow (ogf) (with a proxy compressor power term in the objective) + +"entry point into running the new ogf problem" +function run_nw_ogf(file, model_type, optimizer; kwargs...) + return run_model( + file, + model_type, + optimizer, + build_new_ogf; + solution_processors = [ + sol_psqr_to_p!, + sol_compressor_p_to_r!, + sol_regulator_p_to_r!, + ], + kwargs..., + ) +end + + +"" +function run_soc_new_ogf(file, optimizer; kwargs...) + return run_new_ogf(file, CRDWPGasModel, optimizer; kwargs...) +end + + +"" +function run_dwp_new_ogf(file, optimizer; kwargs...) + return run_new_ogf(file, DWPGasModel, optimizer; kwargs...) +end + + +"construct the new ogf problem" +function build_new_ogf(gm::AbstractGasModel) + bounded_compressors = Dict( + x for x in ref(gm, :compressor) if + _calc_is_compressor_energy_bounded( + get_specific_heat_capacity_ratio(gm.data), + get_gas_specific_gravity(gm.data), + get_temperature(gm.data), + x.second + ) + ) + + variable_pressure(gm) + variable_pressure_sqr(gm) + variable_flow(gm) + variable_on_off_operation(gm) + variable_load_mass_flow(gm) + variable_production_mass_flow(gm) + variable_transfer_mass_flow(gm) + variable_compressor_ratio_sqr(gm) + variable_storage(gm) + variable_form_specific(gm) + + objective_min_new_economic_costs(gm) + + for (i, junction) in ref(gm, :junction) + constraint_mass_flow_balance(gm, i) + + if (junction["junction_type"] == 1) + constraint_pressure(gm, i) + end + end + + for i in ids(gm, :pipe) + constraint_pipe_pressure(gm, i) + constraint_pipe_mass_flow(gm, i) + constraint_pipe_weymouth(gm, i) + end + + for i in ids(gm, :resistor) + constraint_resistor_pressure(gm, i) + constraint_resistor_mass_flow(gm,i) + constraint_resistor_darcy_weisbach(gm,i) + end + + for i in ids(gm, :loss_resistor) + constraint_loss_resistor_pressure(gm, i) + constraint_loss_resistor_mass_flow(gm, i) + end + + for i in ids(gm, :short_pipe) + constraint_short_pipe_pressure(gm, i) + constraint_short_pipe_mass_flow(gm, i) + end + + for i in ids(gm, :compressor) + constraint_compressor_ratios(gm, i) + constraint_compressor_mass_flow(gm, i) + constraint_compressor_ratio_value(gm, i) + end + + for i in keys(bounded_compressors) + constraint_compressor_energy(gm, i) + end + + for i in ids(gm, :valve) + constraint_on_off_valve_mass_flow(gm, i) + constraint_on_off_valve_pressure(gm, i) + end + + for i in ids(gm, :regulator) + constraint_on_off_regulator_mass_flow(gm, i) + constraint_on_off_regulator_pressure(gm, i) + end + +end diff --git a/test/ogf.jl b/test/ogf.jl index 68b25e47..22a6d479 100644 --- a/test/ogf.jl +++ b/test/ogf.jl @@ -10,6 +10,17 @@ @test isapprox(result["solution"]["receipt"]["1"]["fg"], 123.68219958067358; atol = 1e-2) end + @testset "case 6 ogf weymouth lin rel" begin + @info "Testing OGF Linear Relaxation of Pipe Weymouth Physics" + data = GasModels.parse_file("../test/data/matgas/case-6-no-power-limits.m") + result = run_ogf(data, LRWPGasModel, nlp_solver) + @test result["termination_status"] in [LOCALLY_SOLVED, ALMOST_LOCALLY_SOLVED, OPTIMAL, :Suboptimal] + @test isapprox(result["objective"], -260.001; atol = 1e-2) + GasModels.make_si_units!(result["solution"]) + @test isapprox(result["solution"]["receipt"]["1"]["fg"], 130.00040358725565; atol = 1e-2) + end + + @testset "case 6 wp ogf binding energy constraint" begin @info "Testing OGF Binding Energy Cosntraint" data = GasModels.parse_file("../test/data/matgas/case-6.m") From 72d1a777b5a1eedcc8c6e09cf360fa3a65647027 Mon Sep 17 00:00:00 2001 From: Sai Krishna Kanth Hari Date: Tue, 28 Nov 2023 12:32:45 -0600 Subject: [PATCH 3/7] New OGF: Proxy for minimizing compressor power Replacing the nonlinear compressor power objective with a linear proxy, which is the absolute value of the pressure difference at the compressors endpoints --- src/GasModels.jl | 1 + src/core/constraint.jl | 9 ++++++ src/core/constraint_template.jl | 7 +++++ src/core/objective.jl | 4 +-- src/core/variable.jl | 18 ++++++++++++ src/prob/{ogf_new.jl => new_ogf.jl} | 4 ++- test/new_ogf.jl | 44 +++++++++++++++++++++++++++++ test/runtests.jl | 4 ++- 8 files changed, 87 insertions(+), 4 deletions(-) rename src/prob/{ogf_new.jl => new_ogf.jl} (94%) create mode 100644 test/new_ogf.jl diff --git a/src/GasModels.jl b/src/GasModels.jl index 621a9c4c..0591b60e 100644 --- a/src/GasModels.jl +++ b/src/GasModels.jl @@ -83,6 +83,7 @@ module GasModels include("prob/ls.jl") include("prob/nels.jl") include("prob/ogf.jl") + include("prob/new_ogf.jl") include("prob/transient_ogf.jl") include("prob/transient_ogf_archived_storage.jl") diff --git a/src/core/constraint.jl b/src/core/constraint.jl index be1b4760..afb438f4 100644 --- a/src/core/constraint.jl +++ b/src/core/constraint.jl @@ -211,6 +211,15 @@ end function constraint_compressor_mass_flow_ne(gm::AbstractGasModel, n::Int, k, f_min, f_max) end +"Constraint: constraints to serve as a proxy for minimizing the compressor power" +function constraint_compressor_minpower_proxy(gm::AbstractGasModel, n::Int, k, i, j) + pi = var(gm, n, :psqr, i) + pj = var(gm, n, :psqr, j) + mpp = var(gm, n, :min_power_proxy, k) + + _add_constraint!(gm, n, :compressor_minpower_proxy1, k, JuMP.@constraint(gm.model, mpp >= pi - pj )) + _add_constraint!(gm, n, :compressor_minpower_proxy2, k, JuMP.@constraint(gm.model, mpp >= pj - pi )) +end ################################################################################################# # Constraints associated with control valves diff --git a/src/core/constraint_template.jl b/src/core/constraint_template.jl index 71d86e9e..681fb018 100644 --- a/src/core/constraint_template.jl +++ b/src/core/constraint_template.jl @@ -444,7 +444,14 @@ function constraint_compressor_energy(gm::AbstractGasModel, k; n::Int = nw_id_de constraint_compressor_energy(gm, n, k, power_max, m, work, flow_max, ratio_max) end +"Template: Constraints to support the proxy for minimizing the compressor power" +function constraint_compressor_minpower_proxy(gm::AbstractGasModel, k; n::Int = nw_id_default) + compressor = ref(gm, n, :compressor, k) + i = compressor["fr_junction"] + j = compressor["to_junction"] + constraint_compressor_minpower_proxy(gm, n, k, i, j) +end ################################################################################################# # Templates for control valves ################################################################################################# diff --git a/src/core/objective.jl b/src/core/objective.jl index bcabfe99..d2b9d480 100644 --- a/src/core/objective.jl +++ b/src/core/objective.jl @@ -206,7 +206,7 @@ end "function for minimizing new economic costs: ``\\min \\sum_{j \\in {\\cal D}} \\kappa_j \\boldsymbol{d}_j - \\sum_{j \\in {\\cal T}} \\kappa_j \\boldsymbol{\\tau}_j - \\sum_{j \\in {\\cal R}} \\kappa_j \\boldsymbol{r}_j - \\sum_{k \\in {\\cal C}} \\boldsymbol{p^2}_{jk} - \\boldsymbol{p^2}_{ik} ``" function objective_min_new_economic_costs(gm::AbstractGasModel, nws = [nw_id_default]) - r = Dict(n => var(gm, n, :rsqr) for n in nws) + mpp = Dict(n => var(gm, n, :min_power_proxy) for n in nws) f = Dict(n => var(gm, n, :f_compressor) for n in nws) p2 = Dict(n => var(gm, n, :psqr) for n in nws) fl = Dict(n => var(gm, n, :fl) for n in nws) @@ -260,7 +260,7 @@ function objective_min_new_economic_costs(gm::AbstractGasModel, nws = [nw_id_def sum(-transfer_prices[n][i] * ft[n][i] for i in transfer_set[n]) + economic_weighting * sum(prod_prices[n][i] * fg[n][i] for i in prod_set[n]) + (1.0 - economic_weighting) * - sum( (p2[n][compressor["to_junction"]] - p2[n][compressor["fr_junction"]]) for (i, compressor) in ref(gm, n, :compressor)) + sum(mpp[n][i] for (i, compressor) in ref(gm, n, :compressor)) for n in nws )) return JuMP.@objective(gm.model, Min, z) diff --git a/src/core/variable.jl b/src/core/variable.jl index 75ff3f8b..7b2875c3 100644 --- a/src/core/variable.jl +++ b/src/core/variable.jl @@ -685,6 +685,24 @@ function variable_compressor_ratio_sqr_ne(gm::AbstractGasModel, nw::Int=nw_id_de report && sol_component_value(gm, nw, :ne_compressor, :rsqr_ne, keys(compressors), rsqr) end +"Variable Set: variables associated with proxy for minimizing compression power" +function variable_compressor_minpower_proxy(gm::AbstractGasModel, nw::Int=nw_id_default; bounded::Bool=true, report::Bool=true, compressors = ref(gm, nw, :compressor)) + mpp = var(gm, nw)[:min_power_proxy] = JuMP.@variable(gm.model, + [i in keys(compressors)], + base_name="$(nw)_mpp", + start=comp_start_value(ref(gm, nw, :compressor), i, "min_power_proxy_start", 0.0) + ) + + if bounded + for (i, compressor) in compressors + JuMP.set_lower_bound(mpp[i], 0) + end + end + + report && sol_component_value(gm, nw, :compressor, :min_power_proxy, keys(compressors), mpp) +end + + "Support function for getting a one off y direction variable" function get_compressor_y(gm::AbstractGasModel, n::Int, k) if !haskey(var(gm, n),:y_compressor) diff --git a/src/prob/ogf_new.jl b/src/prob/new_ogf.jl similarity index 94% rename from src/prob/ogf_new.jl rename to src/prob/new_ogf.jl index c0c20d1d..73bf1224 100644 --- a/src/prob/ogf_new.jl +++ b/src/prob/new_ogf.jl @@ -1,7 +1,7 @@ # Definitions for running the new optimal gas flow (ogf) (with a proxy compressor power term in the objective) "entry point into running the new ogf problem" -function run_nw_ogf(file, model_type, optimizer; kwargs...) +function run_new_ogf(file, model_type, optimizer; kwargs...) return run_model( file, model_type, @@ -49,6 +49,7 @@ function build_new_ogf(gm::AbstractGasModel) variable_production_mass_flow(gm) variable_transfer_mass_flow(gm) variable_compressor_ratio_sqr(gm) + variable_compressor_minpower_proxy(gm) variable_storage(gm) variable_form_specific(gm) @@ -88,6 +89,7 @@ function build_new_ogf(gm::AbstractGasModel) constraint_compressor_ratios(gm, i) constraint_compressor_mass_flow(gm, i) constraint_compressor_ratio_value(gm, i) + constraint_compressor_minpower_proxy(gm, i) end for i in keys(bounded_compressors) diff --git a/test/new_ogf.jl b/test/new_ogf.jl new file mode 100644 index 00000000..dd1c49d9 --- /dev/null +++ b/test/new_ogf.jl @@ -0,0 +1,44 @@ +@testset "test new ogf" begin + @testset "test wp new ogf" begin + @testset "case 6 new ogf" begin + @info "Testing NEW OGF" + data = GasModels.parse_file("../test/data/matgas/case-6-no-power-limits.m") + result = run_new_ogf(data, WPGasModel, nlp_solver) + @test result["termination_status"] in [LOCALLY_SOLVED, ALMOST_LOCALLY_SOLVED, OPTIMAL, :Suboptimal] + @test isapprox(result["objective"], -253.683; atol = 1e-2) + GasModels.make_si_units!(result["solution"]) + @test isapprox(result["solution"]["receipt"]["1"]["fg"], 123.6822; atol = 1e-2) + end + + @testset "case 6 new ogf weymouth lin rel" begin + @info "Testing NEW OGF Linear Relaxation of Pipe Weymouth Physics" + data = GasModels.parse_file("../test/data/matgas/case-6-no-power-limits.m") + result = run_new_ogf(data, LRWPGasModel, lp_solver) + @test result["termination_status"] in [LOCALLY_SOLVED, ALMOST_LOCALLY_SOLVED, OPTIMAL, :Suboptimal] + @test isapprox(result["objective"], -260.00; atol = 1e-2) + GasModels.make_si_units!(result["solution"]) + @test isapprox(result["solution"]["receipt"]["1"]["fg"], 130.00; atol = 1e-2) + end + + + @testset "case 6 wp new ogf binding energy constraint" begin + @info "Testing NEW OGF Binding Energy Cosntraint" + data = GasModels.parse_file("../test/data/matgas/case-6.m") + result = run_new_ogf(data, WPGasModel, nlp_solver) + @test result["termination_status"] in [LOCALLY_SOLVED, ALMOST_LOCALLY_SOLVED, OPTIMAL, :Suboptimal] + @test isapprox(result["objective"], -237.718; atol = 1e-2) + GasModels.make_si_units!(result["solution"]) + @test isapprox(result["solution"]["receipt"]["1"]["fg"], 97.33; atol = 1e-2) + end + + @testset "case 6 wp new ogf elevation constraint" begin + @info "Testing NEW OGF Elevation Cosntraint" + data = GasModels.parse_file("../test/data/matgas/case-6-elevation.m") + result = run_new_ogf(data, WPGasModel, nlp_solver) + @test result["termination_status"] in [LOCALLY_SOLVED, ALMOST_LOCALLY_SOLVED, OPTIMAL, :Suboptimal] + @test isapprox(result["objective"], -245.251; atol = 1e-2) + GasModels.make_si_units!(result["solution"]) + @test isapprox(result["solution"]["receipt"]["1"]["fg"], 102.5677; atol = 1e-2) + end + end +end diff --git a/test/runtests.jl b/test/runtests.jl index 39f226fc..3bd28ac8 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -17,7 +17,7 @@ import Ipopt import HiGHS import Juniper -using MathOptInterface +using MathOptInterface const MOI = MathOptInterface @@ -52,6 +52,8 @@ include("common.jl") include("ogf.jl") + include("new_ogf.jl") + include("ls.jl") include("nels.jl") From 7a35b047302751d568c3d758268afc82779e867b Mon Sep 17 00:00:00 2001 From: Sai Krishna Kanth Hari Date: Tue, 28 Nov 2023 12:34:24 -0600 Subject: [PATCH 4/7] Description update: New objective Updating the description of the new objective that uses a proxy for the compressor power --- src/core/objective.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/core/objective.jl b/src/core/objective.jl index d2b9d480..4eb41d12 100644 --- a/src/core/objective.jl +++ b/src/core/objective.jl @@ -204,7 +204,7 @@ end "function for minimizing new economic costs: ``\\min \\sum_{j \\in {\\cal D}} \\kappa_j \\boldsymbol{d}_j - \\sum_{j \\in {\\cal T}} \\kappa_j \\boldsymbol{\\tau}_j - \\sum_{j \\in {\\cal R}} \\kappa_j \\boldsymbol{r}_j - - \\sum_{k \\in {\\cal C}} \\boldsymbol{p^2}_{jk} - \\boldsymbol{p^2}_{ik} ``" + \\sum_{k \\in {\\cal C}} |\\boldsymbol{p^2}_{jk} - \\boldsymbol{p^2}_{ik}| ``" function objective_min_new_economic_costs(gm::AbstractGasModel, nws = [nw_id_default]) mpp = Dict(n => var(gm, n, :min_power_proxy) for n in nws) f = Dict(n => var(gm, n, :f_compressor) for n in nws) From 39a5cb79d94fbe70e2b5e022079b1f4397ae0d88 Mon Sep 17 00:00:00 2001 From: Sai Krishna Kanth Hari Date: Wed, 13 Dec 2023 14:50:25 -0600 Subject: [PATCH 5/7] Pipe approximation constraint and new problem Adding a linear approximation for pipe head loss. Replacing the problem ogf_compr_power_proxy by ogf_comp_power_and_pipe_proxy to help the simulator --- src/GasModels.jl | 2 +- src/core/constraint.jl | 16 ++++++++++ src/core/constraint_template.jl | 22 +++++++++++++ ...xy.jl => ogf_comp_power_and_pipe_proxy.jl} | 13 ++++---- ...xy.jl => ogf_comp_power_and_pipe_proxy.jl} | 32 +++++++++---------- test/runtests.jl | 2 +- 6 files changed, 63 insertions(+), 24 deletions(-) rename src/prob/{ogf_comp_power_proxy.jl => ogf_comp_power_and_pipe_proxy.jl} (85%) rename test/{ogf_comp_power_proxy.jl => ogf_comp_power_and_pipe_proxy.jl} (60%) diff --git a/src/GasModels.jl b/src/GasModels.jl index 070a6501..4bb788ec 100644 --- a/src/GasModels.jl +++ b/src/GasModels.jl @@ -83,7 +83,7 @@ module GasModels include("prob/ls.jl") include("prob/nels.jl") include("prob/ogf.jl") - include("prob/ogf_comp_power_proxy.jl") + include("prob/ogf_comp_power_and_pipe_proxy.jl") include("prob/transient_ogf.jl") include("prob/transient_ogf_archived_storage.jl") diff --git a/src/core/constraint.jl b/src/core/constraint.jl index afb438f4..4faaaf18 100644 --- a/src/core/constraint.jl +++ b/src/core/constraint.jl @@ -145,6 +145,22 @@ end ################################################################################################# # Constraints associated with pipes ################################################################################################# +"Constraint: Linear approximation for Weymouth equation" +function constraint_pipe_weymouth_linear_approx(gm::AbstractGasModel, n::Int, k, i, j, f_min, f_max, w, pd_min, pd_max) + pi = var(gm, n, :psqr, i) + pj = var(gm, n, :psqr, j) + f = var(gm, n, :f_pipe, k) + + slope = max(abs(f_max),abs(f_min)) + if w == 0.0 + _add_constraint!(gm, n, :weymouth1, k, JuMP.@constraint(gm.model, (pi - pj) == 0.0)) + elseif f_min == f_max + _add_constraint!(gm, n, :weymouth1, k, JuMP.@constraint(gm.model, w * (pi - pj) == f_min)) + else + _add_constraint!(gm, n, :weymouth1, k, JuMP.@constraint(gm.model, w * (pi - pj) == slope*f)) + end + +end "Constraint: on/off constraints on flow across pipes for expansion pipes" function constraint_pipe_ne(gm::AbstractGasModel, n::Int, k, w, f_min, f_max) diff --git a/src/core/constraint_template.jl b/src/core/constraint_template.jl index 681fb018..554e31df 100644 --- a/src/core/constraint_template.jl +++ b/src/core/constraint_template.jl @@ -146,6 +146,28 @@ function constraint_pipe_weymouth(gm::AbstractGasModel, k; n::Int = nw_id_defaul end end +"Template: Weymouth equation for defining the relationship between pressure drop and flow across a pipe" +function constraint_pipe_weymouth_linear_approx(gm::AbstractGasModel, k; n::Int = nw_id_default) + pipe = ref(gm, n, :pipe, k) + i = pipe["fr_junction"] + j = pipe["to_junction"] + pd_min, pd_max = _calc_pipe_pd_bounds_sqr(pipe, ref(gm, n, :junction, i), ref(gm, n, :junction, j)) + f_min = pipe["flow_min"] + f_max = pipe["flow_max"] + theta = pipe["theta"] + D = pipe["diameter"] + + if(D!=0.0) + if(rad2deg(abs(theta)) <= 5) + w = _calc_pipe_resistance(pipe, gm.ref[:it][gm_it_sym][:base_length], gm.ref[:it][gm_it_sym][:base_pressure], gm.ref[:it][gm_it_sym][:base_flow], gm.ref[:it][gm_it_sym][:sound_speed]) + constraint_pipe_weymouth_linear_approx(gm, n, k, i, j, f_min, f_max, w, pd_min, pd_max) + else + r_1,r_2 = _calc_inclined_pipe_resistance(pipe,gm.ref[:it][gm_it_sym][:base_length], gm.ref[:it][gm_it_sym][:base_pressure], gm.ref[:it][gm_it_sym][:base_flow], gm.ref[:it][gm_it_sym][:sound_speed]) + constraint_inclined_pipe_pressure_drop(gm, n, k, i, j, r_1, r_2) + end + end +end + "Template: Constraint associatd with turning off flow depending on the status of expansion pipes" function constraint_pipe_ne(gm::AbstractGasModel, k; n::Int = nw_id_default) diff --git a/src/prob/ogf_comp_power_proxy.jl b/src/prob/ogf_comp_power_and_pipe_proxy.jl similarity index 85% rename from src/prob/ogf_comp_power_proxy.jl rename to src/prob/ogf_comp_power_and_pipe_proxy.jl index e8db5de8..2cbea805 100644 --- a/src/prob/ogf_comp_power_proxy.jl +++ b/src/prob/ogf_comp_power_and_pipe_proxy.jl @@ -1,12 +1,12 @@ # Definitions for running the new optimal gas flow (ogf) (with a proxy compressor power term in the objective) "entry point into running the new ogf problem" -function run_ogf_comp_power_proxy(file, model_type, optimizer; kwargs...) +function run_ogf_comp_power_and_pipe_proxy(file, model_type, optimizer; kwargs...) return run_model( file, model_type, optimizer, - build_ogf_comp_power_proxy; + build_ogf_comp_power_and_pipe_proxy; solution_processors = [ sol_psqr_to_p!, sol_compressor_p_to_r!, @@ -19,18 +19,18 @@ end "" function run_soc_new_ogf(file, optimizer; kwargs...) - return run_ogf_comp_power_proxy(file, CRDWPGasModel, optimizer; kwargs...) + return run_ogf_comp_power_and_pipe_proxy(file, CRDWPGasModel, optimizer; kwargs...) end "" function run_dwp_new_ogf(file, optimizer; kwargs...) - return run_ogf_comp_power_proxy(file, DWPGasModel, optimizer; kwargs...) + return run_ogf_comp_power_and_pipe_proxy(file, DWPGasModel, optimizer; kwargs...) end "construct the new ogf problem" -function build_ogf_comp_power_proxy(gm::AbstractGasModel) +function build_ogf_comp_power_and_pipe_proxy(gm::AbstractGasModel) bounded_compressors = Dict( x for x in ref(gm, :compressor) if _calc_is_compressor_energy_bounded( @@ -66,7 +66,8 @@ function build_ogf_comp_power_proxy(gm::AbstractGasModel) for i in ids(gm, :pipe) constraint_pipe_pressure(gm, i) constraint_pipe_mass_flow(gm, i) - constraint_pipe_weymouth(gm, i) + # constraint_pipe_weymouth(gm, i) + constraint_pipe_weymouth_linear_approx(gm, i) end for i in ids(gm, :resistor) diff --git a/test/ogf_comp_power_proxy.jl b/test/ogf_comp_power_and_pipe_proxy.jl similarity index 60% rename from test/ogf_comp_power_proxy.jl rename to test/ogf_comp_power_and_pipe_proxy.jl index 1f109a1b..98c63a95 100644 --- a/test/ogf_comp_power_proxy.jl +++ b/test/ogf_comp_power_and_pipe_proxy.jl @@ -1,44 +1,44 @@ @testset "test new ogf" begin @testset "test wp new ogf" begin @testset "case 6 new ogf" begin - @info "Testing OGF with compressor power proxy (linear)" + @info "Testing OGF with compressor power and pipe weymouth proxy (linear)" data = GasModels.parse_file("../test/data/matgas/case-6-no-power-limits.m") - result = run_ogf_comp_power_proxy(data, WPGasModel, nlp_solver) + result = run_ogf_comp_power_and_pipe_proxy(data, WPGasModel, nlp_solver) @test result["termination_status"] in [LOCALLY_SOLVED, ALMOST_LOCALLY_SOLVED, OPTIMAL, :Suboptimal] - @test isapprox(result["objective"], -253.683; atol = 1e-2) + @test isapprox(result["objective"], -230.534; atol = 1e-2) GasModels.make_si_units!(result["solution"]) - @test isapprox(result["solution"]["receipt"]["1"]["fg"], 123.6822; atol = 1e-2) + @test isapprox(result["solution"]["receipt"]["1"]["fg"], 100.53; atol = 1e-2) end @testset "case 6 new ogf weymouth lin rel" begin - @info "Testing OGF with compressor power proxy Linear Relaxation of Pipe Weymouth Physics" + @info "Testing OGF with compressor power proxy and pipe weymouth proxy" data = GasModels.parse_file("../test/data/matgas/case-6-no-power-limits.m") - result = run_ogf_comp_power_proxy(data, LRWPGasModel, lp_solver) + result = run_ogf_comp_power_and_pipe_proxy(data, LRWPGasModel, lp_solver) @test result["termination_status"] in [LOCALLY_SOLVED, ALMOST_LOCALLY_SOLVED, OPTIMAL, :Suboptimal] - @test isapprox(result["objective"], -260.00; atol = 1e-2) + @test isapprox(result["objective"], -230.534; atol = 1e-2) GasModels.make_si_units!(result["solution"]) - @test isapprox(result["solution"]["receipt"]["1"]["fg"], 130.00; atol = 1e-2) + @test isapprox(result["solution"]["receipt"]["1"]["fg"], 100.53; atol = 1e-2) end @testset "case 6 wp new ogf binding energy constraint" begin - @info "Testing OGF with compressor power proxy Binding Energy Cosntraint" + @info "Testing OGF with compressor power and pipe weymouth proxy Binding Energy Cosntraint" data = GasModels.parse_file("../test/data/matgas/case-6.m") - result = run_ogf_comp_power_proxy(data, WPGasModel, nlp_solver) + result = run_ogf_comp_power_and_pipe_proxy(data, WPGasModel, nlp_solver) @test result["termination_status"] in [LOCALLY_SOLVED, ALMOST_LOCALLY_SOLVED, OPTIMAL, :Suboptimal] - @test isapprox(result["objective"], -237.718; atol = 1e-2) + @test isapprox(result["objective"], -187.253; atol = 1e-2) GasModels.make_si_units!(result["solution"]) - @test isapprox(result["solution"]["receipt"]["1"]["fg"], 97.33; atol = 1e-2) + @test isapprox(result["solution"]["receipt"]["1"]["fg"], 67.01; atol = 1e-2) end @testset "case 6 wp new ogf elevation constraint" begin - @info "Testing OGF with compressor power proxy Elevation Cosntraint" + @info "Testing OGF with compressor power and pipe weymouth proxy Elevation Cosntraint" data = GasModels.parse_file("../test/data/matgas/case-6-elevation.m") - result = run_ogf_comp_power_proxy(data, WPGasModel, nlp_solver) + result = run_ogf_comp_power_and_pipe_proxy(data, WPGasModel, nlp_solver) @test result["termination_status"] in [LOCALLY_SOLVED, ALMOST_LOCALLY_SOLVED, OPTIMAL, :Suboptimal] - @test isapprox(result["objective"], -245.251; atol = 1e-2) + @test isapprox(result["objective"], -219.461; atol = 1e-2) GasModels.make_si_units!(result["solution"]) - @test isapprox(result["solution"]["receipt"]["1"]["fg"], 102.5677; atol = 1e-2) + @test isapprox(result["solution"]["receipt"]["1"]["fg"], 86.33; atol = 1e-2) end end end diff --git a/test/runtests.jl b/test/runtests.jl index 5b3ba822..bce31467 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -52,7 +52,7 @@ include("common.jl") include("ogf.jl") - include("ogf_comp_power_proxy.jl") + include("ogf_comp_power_and_pipe_proxy.jl") include("ls.jl") From 133035e2b8a9ee0e46d442de32127fcc1550c9be Mon Sep 17 00:00:00 2001 From: Sai Krishna Kanth Hari Date: Wed, 13 Dec 2023 14:54:45 -0600 Subject: [PATCH 6/7] removing old files --- src/prob/new_ogf.jl | 109 -------------------------------------------- test/new_ogf.jl | 44 ------------------ 2 files changed, 153 deletions(-) delete mode 100644 src/prob/new_ogf.jl delete mode 100644 test/new_ogf.jl diff --git a/src/prob/new_ogf.jl b/src/prob/new_ogf.jl deleted file mode 100644 index 73bf1224..00000000 --- a/src/prob/new_ogf.jl +++ /dev/null @@ -1,109 +0,0 @@ -# Definitions for running the new optimal gas flow (ogf) (with a proxy compressor power term in the objective) - -"entry point into running the new ogf problem" -function run_new_ogf(file, model_type, optimizer; kwargs...) - return run_model( - file, - model_type, - optimizer, - build_new_ogf; - solution_processors = [ - sol_psqr_to_p!, - sol_compressor_p_to_r!, - sol_regulator_p_to_r!, - ], - kwargs..., - ) -end - - -"" -function run_soc_new_ogf(file, optimizer; kwargs...) - return run_new_ogf(file, CRDWPGasModel, optimizer; kwargs...) -end - - -"" -function run_dwp_new_ogf(file, optimizer; kwargs...) - return run_new_ogf(file, DWPGasModel, optimizer; kwargs...) -end - - -"construct the new ogf problem" -function build_new_ogf(gm::AbstractGasModel) - bounded_compressors = Dict( - x for x in ref(gm, :compressor) if - _calc_is_compressor_energy_bounded( - get_specific_heat_capacity_ratio(gm.data), - get_gas_specific_gravity(gm.data), - get_temperature(gm.data), - x.second - ) - ) - - variable_pressure(gm) - variable_pressure_sqr(gm) - variable_flow(gm) - variable_on_off_operation(gm) - variable_load_mass_flow(gm) - variable_production_mass_flow(gm) - variable_transfer_mass_flow(gm) - variable_compressor_ratio_sqr(gm) - variable_compressor_minpower_proxy(gm) - variable_storage(gm) - variable_form_specific(gm) - - objective_min_new_economic_costs(gm) - - for (i, junction) in ref(gm, :junction) - constraint_mass_flow_balance(gm, i) - - if (junction["junction_type"] == 1) - constraint_pressure(gm, i) - end - end - - for i in ids(gm, :pipe) - constraint_pipe_pressure(gm, i) - constraint_pipe_mass_flow(gm, i) - constraint_pipe_weymouth(gm, i) - end - - for i in ids(gm, :resistor) - constraint_resistor_pressure(gm, i) - constraint_resistor_mass_flow(gm,i) - constraint_resistor_darcy_weisbach(gm,i) - end - - for i in ids(gm, :loss_resistor) - constraint_loss_resistor_pressure(gm, i) - constraint_loss_resistor_mass_flow(gm, i) - end - - for i in ids(gm, :short_pipe) - constraint_short_pipe_pressure(gm, i) - constraint_short_pipe_mass_flow(gm, i) - end - - for i in ids(gm, :compressor) - constraint_compressor_ratios(gm, i) - constraint_compressor_mass_flow(gm, i) - constraint_compressor_ratio_value(gm, i) - constraint_compressor_minpower_proxy(gm, i) - end - - for i in keys(bounded_compressors) - constraint_compressor_energy(gm, i) - end - - for i in ids(gm, :valve) - constraint_on_off_valve_mass_flow(gm, i) - constraint_on_off_valve_pressure(gm, i) - end - - for i in ids(gm, :regulator) - constraint_on_off_regulator_mass_flow(gm, i) - constraint_on_off_regulator_pressure(gm, i) - end - -end diff --git a/test/new_ogf.jl b/test/new_ogf.jl deleted file mode 100644 index dd1c49d9..00000000 --- a/test/new_ogf.jl +++ /dev/null @@ -1,44 +0,0 @@ -@testset "test new ogf" begin - @testset "test wp new ogf" begin - @testset "case 6 new ogf" begin - @info "Testing NEW OGF" - data = GasModels.parse_file("../test/data/matgas/case-6-no-power-limits.m") - result = run_new_ogf(data, WPGasModel, nlp_solver) - @test result["termination_status"] in [LOCALLY_SOLVED, ALMOST_LOCALLY_SOLVED, OPTIMAL, :Suboptimal] - @test isapprox(result["objective"], -253.683; atol = 1e-2) - GasModels.make_si_units!(result["solution"]) - @test isapprox(result["solution"]["receipt"]["1"]["fg"], 123.6822; atol = 1e-2) - end - - @testset "case 6 new ogf weymouth lin rel" begin - @info "Testing NEW OGF Linear Relaxation of Pipe Weymouth Physics" - data = GasModels.parse_file("../test/data/matgas/case-6-no-power-limits.m") - result = run_new_ogf(data, LRWPGasModel, lp_solver) - @test result["termination_status"] in [LOCALLY_SOLVED, ALMOST_LOCALLY_SOLVED, OPTIMAL, :Suboptimal] - @test isapprox(result["objective"], -260.00; atol = 1e-2) - GasModels.make_si_units!(result["solution"]) - @test isapprox(result["solution"]["receipt"]["1"]["fg"], 130.00; atol = 1e-2) - end - - - @testset "case 6 wp new ogf binding energy constraint" begin - @info "Testing NEW OGF Binding Energy Cosntraint" - data = GasModels.parse_file("../test/data/matgas/case-6.m") - result = run_new_ogf(data, WPGasModel, nlp_solver) - @test result["termination_status"] in [LOCALLY_SOLVED, ALMOST_LOCALLY_SOLVED, OPTIMAL, :Suboptimal] - @test isapprox(result["objective"], -237.718; atol = 1e-2) - GasModels.make_si_units!(result["solution"]) - @test isapprox(result["solution"]["receipt"]["1"]["fg"], 97.33; atol = 1e-2) - end - - @testset "case 6 wp new ogf elevation constraint" begin - @info "Testing NEW OGF Elevation Cosntraint" - data = GasModels.parse_file("../test/data/matgas/case-6-elevation.m") - result = run_new_ogf(data, WPGasModel, nlp_solver) - @test result["termination_status"] in [LOCALLY_SOLVED, ALMOST_LOCALLY_SOLVED, OPTIMAL, :Suboptimal] - @test isapprox(result["objective"], -245.251; atol = 1e-2) - GasModels.make_si_units!(result["solution"]) - @test isapprox(result["solution"]["receipt"]["1"]["fg"], 102.5677; atol = 1e-2) - end - end -end From 45d10393d51c82fb137a7ba7ba4f708b81226f16 Mon Sep 17 00:00:00 2001 From: ftuser Date: Thu, 14 Dec 2023 13:13:11 -0700 Subject: [PATCH 7/7] adds linear approx --- CHANGELOG.md | 3 + src/GasModels.jl | 1 + src/prob/ogf_comp_power_and_pipe_proxy.jl | 4 +- src/prob/ogf_comp_power_proxy.jl | 109 ++++++++++++++++++++++ test/ogf_comp_power_proxy.jl | 44 +++++++++ test/runtests.jl | 2 + 6 files changed, 161 insertions(+), 2 deletions(-) create mode 100644 src/prob/ogf_comp_power_proxy.jl create mode 100644 test/ogf_comp_power_proxy.jl diff --git a/CHANGELOG.md b/CHANGELOG.md index f32a2110..c8ca6513 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -8,6 +8,9 @@ - adds linear relaxation formulation for steady state OGF - adds a new `run_ogf_comp_power_proxy` problem formulation - the `LRWPGasModel` now is populated +- adds linear approximation formulation for steady state OGF +- adds a new `run_ogf_comp_power_and_pipe_proxy` problem formulation +- the `LRWPGasModel` now is populated ## v0.9.2 diff --git a/src/GasModels.jl b/src/GasModels.jl index 4bb788ec..d5d754a1 100644 --- a/src/GasModels.jl +++ b/src/GasModels.jl @@ -83,6 +83,7 @@ module GasModels include("prob/ls.jl") include("prob/nels.jl") include("prob/ogf.jl") + include("prob/ogf_comp_power_proxy.jl") include("prob/ogf_comp_power_and_pipe_proxy.jl") include("prob/transient_ogf.jl") include("prob/transient_ogf_archived_storage.jl") diff --git a/src/prob/ogf_comp_power_and_pipe_proxy.jl b/src/prob/ogf_comp_power_and_pipe_proxy.jl index 2cbea805..79c96d64 100644 --- a/src/prob/ogf_comp_power_and_pipe_proxy.jl +++ b/src/prob/ogf_comp_power_and_pipe_proxy.jl @@ -18,13 +18,13 @@ end "" -function run_soc_new_ogf(file, optimizer; kwargs...) +function run_soc_ogf_comp_power_and_pipe_proxy(file, optimizer; kwargs...) return run_ogf_comp_power_and_pipe_proxy(file, CRDWPGasModel, optimizer; kwargs...) end "" -function run_dwp_new_ogf(file, optimizer; kwargs...) +function run_dwp_ogf_comp_power_and_pipe_proxy(file, optimizer; kwargs...) return run_ogf_comp_power_and_pipe_proxy(file, DWPGasModel, optimizer; kwargs...) end diff --git a/src/prob/ogf_comp_power_proxy.jl b/src/prob/ogf_comp_power_proxy.jl new file mode 100644 index 00000000..477fa086 --- /dev/null +++ b/src/prob/ogf_comp_power_proxy.jl @@ -0,0 +1,109 @@ +# Definitions for running the new optimal gas flow (ogf) (with a proxy compressor power term in the objective) + +"entry point into running the new ogf problem" +function run_ogf_comp_power_proxy(file, model_type, optimizer; kwargs...) + return run_model( + file, + model_type, + optimizer, + build_ogf_comp_power_proxy; + solution_processors = [ + sol_psqr_to_p!, + sol_compressor_p_to_r!, + sol_regulator_p_to_r!, + ], + kwargs..., + ) +end + + +"" +function run_soc_ogf_comp_power_proxy(file, optimizer; kwargs...) + return run_ogf_comp_power_proxy(file, CRDWPGasModel, optimizer; kwargs...) +end + + +"" +function run_dwp_ogf_comp_power_proxy(file, optimizer; kwargs...) + return run_ogf_comp_power_proxy(file, DWPGasModel, optimizer; kwargs...) +end + + +"construct the new ogf problem" +function build_ogf_comp_power_proxy(gm::AbstractGasModel) + bounded_compressors = Dict( + x for x in ref(gm, :compressor) if + _calc_is_compressor_energy_bounded( + get_specific_heat_capacity_ratio(gm.data), + get_gas_specific_gravity(gm.data), + get_temperature(gm.data), + x.second + ) + ) + + variable_pressure(gm) + variable_pressure_sqr(gm) + variable_flow(gm) + variable_on_off_operation(gm) + variable_load_mass_flow(gm) + variable_production_mass_flow(gm) + variable_transfer_mass_flow(gm) + variable_compressor_ratio_sqr(gm) + variable_compressor_minpower_proxy(gm) + variable_storage(gm) + variable_form_specific(gm) + + objective_min_proxy_economic_costs(gm) + + for (i, junction) in ref(gm, :junction) + constraint_mass_flow_balance(gm, i) + + if (junction["junction_type"] == 1) + constraint_pressure(gm, i) + end + end + + for i in ids(gm, :pipe) + constraint_pipe_pressure(gm, i) + constraint_pipe_mass_flow(gm, i) + constraint_pipe_weymouth(gm, i) + end + + for i in ids(gm, :resistor) + constraint_resistor_pressure(gm, i) + constraint_resistor_mass_flow(gm,i) + constraint_resistor_darcy_weisbach(gm,i) + end + + for i in ids(gm, :loss_resistor) + constraint_loss_resistor_pressure(gm, i) + constraint_loss_resistor_mass_flow(gm, i) + end + + for i in ids(gm, :short_pipe) + constraint_short_pipe_pressure(gm, i) + constraint_short_pipe_mass_flow(gm, i) + end + + for i in ids(gm, :compressor) + constraint_compressor_ratios(gm, i) + constraint_compressor_mass_flow(gm, i) + constraint_compressor_ratio_value(gm, i) + constraint_compressor_minpower_proxy(gm, i) + end + + for i in keys(bounded_compressors) + constraint_compressor_energy(gm, i) + end + + for i in ids(gm, :valve) + constraint_on_off_valve_mass_flow(gm, i) + constraint_on_off_valve_pressure(gm, i) + end + + for i in ids(gm, :regulator) + constraint_on_off_regulator_mass_flow(gm, i) + constraint_on_off_regulator_pressure(gm, i) + end + +end \ No newline at end of file diff --git a/test/ogf_comp_power_proxy.jl b/test/ogf_comp_power_proxy.jl new file mode 100644 index 00000000..37628ed3 --- /dev/null +++ b/test/ogf_comp_power_proxy.jl @@ -0,0 +1,44 @@ +@testset "test new ogf" begin + @testset "test wp new ogf" begin + @testset "case 6 new ogf" begin + @info "Testing OGF with compressor power proxy (linear)" + data = GasModels.parse_file("../test/data/matgas/case-6-no-power-limits.m") + result = run_ogf_comp_power_proxy(data, WPGasModel, nlp_solver) + @test result["termination_status"] in [LOCALLY_SOLVED, ALMOST_LOCALLY_SOLVED, OPTIMAL, :Suboptimal] + @test isapprox(result["objective"], -253.683; atol = 1e-2) + GasModels.make_si_units!(result["solution"]) + @test isapprox(result["solution"]["receipt"]["1"]["fg"], 123.6822; atol = 1e-2) + end + + @testset "case 6 new ogf weymouth lin rel" begin + @info "Testing OGF with compressor power proxy Linear Relaxation of Pipe Weymouth Physics" + data = GasModels.parse_file("../test/data/matgas/case-6-no-power-limits.m") + result = run_ogf_comp_power_proxy(data, LRWPGasModel, lp_solver) + @test result["termination_status"] in [LOCALLY_SOLVED, ALMOST_LOCALLY_SOLVED, OPTIMAL, :Suboptimal] + @test isapprox(result["objective"], -260.00; atol = 1e-2) + GasModels.make_si_units!(result["solution"]) + @test isapprox(result["solution"]["receipt"]["1"]["fg"], 130.00; atol = 1e-2) + end + + + @testset "case 6 wp new ogf binding energy constraint" begin + @info "Testing OGF with compressor power proxy Binding Energy Cosntraint" + data = GasModels.parse_file("../test/data/matgas/case-6.m") + result = run_ogf_comp_power_proxy(data, WPGasModel, nlp_solver) + @test result["termination_status"] in [LOCALLY_SOLVED, ALMOST_LOCALLY_SOLVED, OPTIMAL, :Suboptimal] + @test isapprox(result["objective"], -237.718; atol = 1e-2) + GasModels.make_si_units!(result["solution"]) + @test isapprox(result["solution"]["receipt"]["1"]["fg"], 97.33; atol = 1e-2) + end + + @testset "case 6 wp new ogf elevation constraint" begin + @info "Testing OGF with compressor power proxy Elevation Cosntraint" + data = GasModels.parse_file("../test/data/matgas/case-6-elevation.m") + result = run_ogf_comp_power_proxy(data, WPGasModel, nlp_solver) + @test result["termination_status"] in [LOCALLY_SOLVED, ALMOST_LOCALLY_SOLVED, OPTIMAL, :Suboptimal] + @test isapprox(result["objective"], -245.251; atol = 1e-2) + GasModels.make_si_units!(result["solution"]) + @test isapprox(result["solution"]["receipt"]["1"]["fg"], 102.5677; atol = 1e-2) + end + end +end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index bce31467..f92f38f7 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -54,6 +54,8 @@ include("common.jl") include("ogf_comp_power_and_pipe_proxy.jl") + include("ogf_comp_power_proxy.jl") + include("ls.jl") include("nels.jl")