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

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

src/GasModels.jl

@@ -12,6 +12,8 @@ module GasModels
 
     using Dates
     using Dierckx
+    using PolyhedralRelaxations
+
 
 
     # Create our module level logger (this will get precompiled)

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

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)

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)

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)

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
 
 

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)

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)

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)
 

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)

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)

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)
 

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

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