Merge pull request #255 from lanl-ansi/linearapprox

Linear Approximation

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 

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")
 

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)

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)

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

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

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

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

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")