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

Author: hskkanth

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

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

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
 #################################################################################################

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)

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)

src/prob/ogf_new.jl renamed to 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)

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

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