Add linear branch flow model and branch flow storage problem type

CHANGELOG.md

@@ -2,7 +2,8 @@ PowerModels.jl Change Log
 =========================
 
 ### Staged
-- nothing
+- Add support for building storage models with branch flow formulations (#676)
+- A linear branch flow formulation was added (#676)
 
 ### v0.15.2
 - Add support for on/off storage in SOCWRPowerModel

README.md

@@ -65,6 +65,7 @@ The primary developer is Carleton Coffrin(@ccoffrin) with support from the follo
 - Miles Lubin (@mlubin) MIT, Julia/JuMP advise
 - Yeesian Ng (@yeesian) MIT, Documenter.jl setup
 - Kaarthik Sundar (@kaarthiksundar) LANL, OBBT utility
+- Per Aaslid (@peraaslid) SINTEF ER, Branch flow storage model and linear branch flow formulation
 
 
 ## Citing PowerModels

src/core/types.jl

@@ -28,6 +28,10 @@ abstract type AbstractConicModel <: AbstractPowerModel end
 "for branch flow models"
 abstract type AbstractBFModel <: AbstractPowerModel end
 
+
+"for variants of branch flow models that target LP solvers"
+abstract type AbstractBFAModel <: AbstractBFModel end
+
 "for variants of branch flow models that target QP or NLP solvers"
 abstract type AbstractBFQPModel <: AbstractBFModel end
 
@@ -389,6 +393,28 @@ Extended as discussed in:
 """
 mutable struct SOCBFPowerModel <: AbstractSOCBFModel @pm_fields end
 
+"""
+Linear approximation of branch flow model.
+
+The implementation builds on the second-order cone relaxation of the branch
+flow model, but neglects the active and reactive loss terms associated with
+the squared current magnitude so the power flow equations become linear.
+Note that flow bounds are still second order cones.
+```
+@article{Baran1989OptimalSystems,
+    title = {{Optimal capacitor placement on radial distribution systems}},
+    year = {1989},
+    journal = {IEEE Transactions on Power Delivery},
+    author = {Baran, Mesut E. and Wu, Felix F.},
+    number = {1},
+    pages = {725--734},
+    volume = {4},
+    doi = {10.1109/61.19265},
+    issn = {19374208}
+}
+```
+"""
+mutable struct BFAPowerModel <: AbstractBFAModel @pm_fields end
 
 ""
 abstract type AbstractSOCBFConicModel <: AbstractBFConicModel end

src/form/bf.jl

@@ -150,3 +150,64 @@ function constraint_voltage_angle_difference(pm::AbstractBFModel, n::Int, f_idx,
                  >= ((ti + tzi*g_fr - tzr*b_fr)*(w_fr/tm^2) - tzi*p_fr - tzr*q_fr)
         )
 end
+
+"""
+Defines linear branch flow model power flow equations
+"""
+function constraint_flow_losses(pm::AbstractBFAModel, n::Int, i, f_bus, t_bus, f_idx, t_idx, r, x, g_sh_fr, g_sh_to, b_sh_fr, b_sh_to, tm)
+    p_fr = var(pm, n, :p, f_idx)
+    q_fr = var(pm, n, :q, f_idx)
+    p_to = var(pm, n, :p, t_idx)
+    q_to = var(pm, n, :q, t_idx)
+    w_fr = var(pm, n, :w, f_bus)
+    w_to = var(pm, n, :w, t_bus)
+
+    JuMP.@constraint(pm.model, p_fr + p_to == g_sh_fr*(w_fr/tm^2) + g_sh_to*w_to)
+    JuMP.@constraint(pm.model, q_fr + q_to == b_sh_fr*(w_fr/tm^2) + b_sh_to*w_to)
+end
+
+
+"""
+Defines voltage drop over a branch, linking from and to side voltage magnitude
+"""
+function constraint_voltage_magnitude_difference(pm::AbstractBFAModel, n::Int, i, f_bus, t_bus, f_idx, t_idx, r, x, g_sh_fr, b_sh_fr, tm)
+    p_fr = var(pm, n, :p, f_idx)
+    q_fr = var(pm, n, :q, f_idx)
+    w_fr = var(pm, n, :w, f_bus)
+    w_to = var(pm, n, :w, t_bus)
+
+    JuMP.@constraint(pm.model, (w_fr/tm^2) - w_to ==  2*(r*p_fr + x*q_fr))
+end
+
+
+""
+function constraint_model_current(pm::AbstractBFAModel, n::Int)
+
+end
+
+""
+function variable_current_magnitude_sqr(pm::AbstractBFAModel; nw::Int=pm.cnw, bounded::Bool=true, report::Bool=true)
+
+end
+
+
+"Neglects the active and reactive loss terms associated with the squared current magnitude."
+function constraint_storage_loss(pm::AbstractBFAModel, n::Int, i, bus, r, x, p_loss, q_loss; conductors=[1])
+    ps = var(pm, n, :ps, i)
+    qs = var(pm, n, :qs, i)
+    sc = var(pm, n, :sc, i)
+    sd = var(pm, n, :sd, i)
+
+
+    JuMP.@constraint(pm.model,
+        sum(ps[c] for c in conductors) + (sd - sc)
+        ==
+        p_loss
+    )
+
+    JuMP.@constraint(pm.model,
+        sum(qs[c] for c in conductors)
+        ==
+        q_loss
+    )
+end

src/form/shared.jl

@@ -15,14 +15,14 @@
 function variable_shunt_factor(pm::AbstractWConvexModels; nw::Int=pm.cnw, relax::Bool=false, report::Bool=true)
     if !relax
         z_shunt = var(pm, nw)[:z_shunt] = JuMP.@variable(pm.model,
-            [i in ids(pm, nw, :shunt)], base_name="$(nw)_z_shunt", 
+            [i in ids(pm, nw, :shunt)], base_name="$(nw)_z_shunt",
             binary = true,
             start = comp_start_value(ref(pm, nw, :shunt, i), "z_shunt_start", 1.0)
         )
     else
         z_shunt = var(pm, nw)[:z_shunt] = JuMP.@variable(pm.model,
-            [i in ids(pm, nw, :shunt)], base_name="$(nw)_z_shunt", 
-            upper_bound = 1, 
+            [i in ids(pm, nw, :shunt)], base_name="$(nw)_z_shunt",
+            upper_bound = 1,
             lower_bound = 0,
             start = comp_start_value(ref(pm, nw, :shunt, i), "z_shunt_start", 1.0)
         )
@@ -39,6 +39,33 @@ function variable_shunt_factor(pm::AbstractWConvexModels; nw::Int=pm.cnw, relax:
 end
 
 
+"do nothing by default but some formulations require this"
+function variable_current_storage(pm::AbstractWConvexModels; nw::Int=pm.cnw, bounded::Bool=true, report::Bool=true)
+    ccms = var(pm, nw)[:ccms] = JuMP.@variable(pm.model,
+        [i in ids(pm, nw, :storage)], base_name="$(nw)_ccms",
+        start = comp_start_value(ref(pm, nw, :storage, i), "ccms_start")
+    )
+
+    if bounded
+        bus = ref(pm, nw, :bus)
+        for (i, storage) in ref(pm, nw, :storage)
+            ub = Inf
+            if haskey(storage, "thermal_rating")
+                sb = bus[storage["storage_bus"]]
+                ub = (storage["thermal_rating"]/sb["vmin"])^2
+            end
+
+            JuMP.set_lower_bound(ccms[i], 0.0)
+            if !isinf(ub)
+                JuMP.set_upper_bound(ccms[i], ub)
+            end
+        end
+    end
+
+    report && sol_component_value(pm, nw, :storage, :ccms, ids(pm, nw, :storage), ccms)
+end
+
+
 "`t[ref_bus] == 0`"
 function constraint_theta_ref(pm::AbstractPolarModels, n::Int, i::Int)
     JuMP.@constraint(pm.model, var(pm, n, :va)[i] == 0)
@@ -343,3 +370,28 @@ function constraint_current_limit(pm::AbstractWModels, n::Int, f_idx, c_rating_a
     JuMP.@constraint(pm.model, p_to^2 + q_to^2 <= w_to*c_rating_a^2)
 end
 
+""
+function constraint_storage_loss(pm::AbstractWConvexModels, n::Int, i, bus, r, x, p_loss, q_loss; conductors=[1])
+    w = var(pm, n, :w, bus)
+    ccms = var(pm, n, :ccms, i)
+    ps = var(pm, n, :ps, i)
+    qs = var(pm, n, :qs, i)
+    sc = var(pm, n, :sc, i)
+    sd = var(pm, n, :sd, i)
+
+    for c in conductors
+        JuMP.@constraint(pm.model, ps[c]^2 + qs[c]^2 <= w[c]*ccms[c])
+    end
+
+    JuMP.@constraint(pm.model,
+        sum(ps[c] for c in conductors) + (sd - sc)
+        ==
+        p_loss + sum(r[c]*ccms[c] for c in conductors)
+    )
+
+    JuMP.@constraint(pm.model,
+        sum(qs[c] for c in conductors)
+        ==
+        q_loss + sum(x[c]*ccms[c] for c in conductors)
+    )
+end

src/form/wr.jl

@@ -336,60 +336,6 @@ function variable_voltage_product_ne(pm::AbstractWRModel; nw::Int=pm.cnw, report
 end
 
 
-"do nothing by default but some formulations require this"
-function variable_current_storage(pm::AbstractWRModel; nw::Int=pm.cnw, bounded::Bool=true, report::Bool=true)
-    ccms = var(pm, nw)[:ccms] = JuMP.@variable(pm.model,
-        [i in ids(pm, nw, :storage)], base_name="$(nw)_ccms",
-        start = comp_start_value(ref(pm, nw, :storage, i), "ccms_start")
-    )
-
-    if bounded
-        bus = ref(pm, nw, :bus)
-        for (i, storage) in ref(pm, nw, :storage)
-            ub = Inf
-            if haskey(storage, "thermal_rating")
-                sb = bus[storage["storage_bus"]]
-                ub = (storage["thermal_rating"]/sb["vmin"])^2
-            end
-
-            JuMP.set_lower_bound(ccms[i], 0.0)
-            if !isinf(ub)
-                JuMP.set_upper_bound(ccms[i], ub)
-            end
-        end
-    end
-
-    report && sol_component_value(pm, nw, :storage, :ccms, ids(pm, nw, :storage), ccms)
-end
-
-""
-function constraint_storage_loss(pm::AbstractWRModel, n::Int, i, bus, r, x, p_loss, q_loss; conductors=[1])
-    w = var(pm, n, :w, bus)
-    ccms = var(pm, n, :ccms, i)
-    ps = var(pm, n, :ps, i)
-    qs = var(pm, n, :qs, i)
-    sc = var(pm, n, :sc, i)
-    sd = var(pm, n, :sd, i)
-
-    for c in conductors
-        JuMP.@constraint(pm.model, ps[c]^2 + qs[c]^2 <= w[c]*ccms[c])
-    end
-
-    JuMP.@constraint(pm.model, 
-        sum(ps[c] for c in conductors) + (sd - sc)
-        ==
-        p_loss + sum(r[c]*ccms[c] for c in conductors)
-    )
-
-    JuMP.@constraint(pm.model, 
-        sum(qs[c] for c in conductors)
-        ==
-        q_loss + sum(x[c]*ccms[c] for c in conductors)
-    )
-end
-
-
-
 ###### QC Relaxations ######
 
 

src/prob/opf_bf.jl

@@ -3,6 +3,11 @@ function run_opf_bf(file, model_type::Type{T}, optimizer; kwargs...) where T <:
     return run_model(file, model_type, optimizer, build_opf_bf; kwargs...)
 end
 
+""
+function run_mn_opf_bf_strg(file, model_type::Type, optimizer; kwargs...)
+    return run_model(file, model_type, optimizer, build_mn_opf_bf_strg; multinetwork=true, kwargs...)
+end
+
 ""
 function build_opf_bf(pm::AbstractPowerModel)
     variable_voltage(pm)
@@ -37,3 +42,61 @@ function build_opf_bf(pm::AbstractPowerModel)
         constraint_dcline(pm, i)
     end
 end
+
+"Build multinetwork branch flow storage OPF"
+function build_mn_opf_bf_strg(pm::AbstractPowerModel)
+    for (n, network) in nws(pm)
+        variable_voltage(pm, nw=n)
+        variable_generation(pm, nw=n)
+        variable_storage_mi(pm, nw=n)
+        variable_branch_flow(pm, nw=n)
+        variable_branch_current(pm, nw=n)
+        variable_dcline_flow(pm, nw=n)
+
+        constraint_model_current(pm, nw=n)
+
+        for i in ids(pm, :ref_buses, nw=n)
+            constraint_theta_ref(pm, i, nw=n)
+        end
+
+        for i in ids(pm, :bus, nw=n)
+            constraint_power_balance(pm, i, nw=n)
+        end
+
+        for i in ids(pm, :storage, nw=n)
+            constraint_storage_complementarity_mi(pm, i, nw=n)
+            constraint_storage_loss(pm, i, nw=n)
+            constraint_storage_thermal_limit(pm, i, nw=n)
+        end
+
+        for i in ids(pm, :branch, nw=n)
+            constraint_flow_losses(pm, i, nw=n)
+            constraint_voltage_magnitude_difference(pm, i, nw=n)
+
+            constraint_voltage_angle_difference(pm, i, nw=n)
+
+            constraint_thermal_limit_from(pm, i, nw=n)
+            constraint_thermal_limit_to(pm, i, nw=n)
+        end
+
+        for i in ids(pm, :dcline, nw=n)
+            constraint_dcline(pm, i, nw=n)
+        end
+    end
+
+    network_ids = sort(collect(nw_ids(pm)))
+
+    n_1 = network_ids[1]
+    for i in ids(pm, :storage, nw=n_1)
+        constraint_storage_state(pm, i, nw=n_1)
+    end
+
+    for n_2 in network_ids[2:end]
+        for i in ids(pm, :storage, nw=n_2)
+            constraint_storage_state(pm, i, n_1, n_2)
+        end
+        n_1 = n_2
+    end
+
+    objective_min_fuel_and_flow_cost(pm)
+end