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 070a6501..d5d754a1 100644 --- a/src/GasModels.jl +++ b/src/GasModels.jl @@ -84,6 +84,7 @@ module GasModels 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_and_pipe_proxy.jl b/src/prob/ogf_comp_power_and_pipe_proxy.jl new file mode 100644 index 00000000..79c96d64 --- /dev/null +++ b/src/prob/ogf_comp_power_and_pipe_proxy.jl @@ -0,0 +1,110 @@ +# 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_and_pipe_proxy(file, model_type, optimizer; kwargs...) + return run_model( + file, + model_type, + optimizer, + build_ogf_comp_power_and_pipe_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_and_pipe_proxy(file, optimizer; kwargs...) + return run_ogf_comp_power_and_pipe_proxy(file, CRDWPGasModel, optimizer; kwargs...) +end + + +"" +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 + + +"construct the new ogf problem" +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( + 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) + constraint_pipe_weymouth_linear_approx(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/src/prob/ogf_comp_power_proxy.jl b/src/prob/ogf_comp_power_proxy.jl index e8db5de8..477fa086 100644 --- a/src/prob/ogf_comp_power_proxy.jl +++ b/src/prob/ogf_comp_power_proxy.jl @@ -18,13 +18,13 @@ end "" -function run_soc_new_ogf(file, optimizer; kwargs...) +function run_soc_ogf_comp_power_proxy(file, optimizer; kwargs...) return run_ogf_comp_power_proxy(file, CRDWPGasModel, optimizer; kwargs...) end "" -function run_dwp_new_ogf(file, optimizer; kwargs...) +function run_dwp_ogf_comp_power_proxy(file, optimizer; kwargs...) return run_ogf_comp_power_proxy(file, DWPGasModel, optimizer; kwargs...) end @@ -106,4 +106,4 @@ function build_ogf_comp_power_proxy(gm::AbstractGasModel) constraint_on_off_regulator_pressure(gm, i) end -end +end \ No newline at end of file diff --git a/test/ogf_comp_power_and_pipe_proxy.jl b/test/ogf_comp_power_and_pipe_proxy.jl new file mode 100644 index 00000000..98c63a95 --- /dev/null +++ b/test/ogf_comp_power_and_pipe_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 and pipe weymouth proxy (linear)" + data = GasModels.parse_file("../test/data/matgas/case-6-no-power-limits.m") + 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"], -230.534; atol = 1e-2) + GasModels.make_si_units!(result["solution"]) + @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 and pipe weymouth proxy" + data = GasModels.parse_file("../test/data/matgas/case-6-no-power-limits.m") + 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"], -230.534; atol = 1e-2) + GasModels.make_si_units!(result["solution"]) + @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 and pipe weymouth proxy Binding Energy Cosntraint" + data = GasModels.parse_file("../test/data/matgas/case-6.m") + 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"], -187.253; atol = 1e-2) + GasModels.make_si_units!(result["solution"]) + @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 and pipe weymouth proxy Elevation Cosntraint" + data = GasModels.parse_file("../test/data/matgas/case-6-elevation.m") + 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"], -219.461; atol = 1e-2) + GasModels.make_si_units!(result["solution"]) + @test isapprox(result["solution"]["receipt"]["1"]["fg"], 86.33; atol = 1e-2) + end + end +end diff --git a/test/ogf_comp_power_proxy.jl b/test/ogf_comp_power_proxy.jl index 1f109a1b..37628ed3 100644 --- a/test/ogf_comp_power_proxy.jl +++ b/test/ogf_comp_power_proxy.jl @@ -41,4 +41,4 @@ @test isapprox(result["solution"]["receipt"]["1"]["fg"], 102.5677; atol = 1e-2) end end -end +end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 5b3ba822..f92f38f7 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -52,6 +52,8 @@ include("common.jl") include("ogf.jl") + include("ogf_comp_power_and_pipe_proxy.jl") + include("ogf_comp_power_proxy.jl") include("ls.jl")