Add blocker placement and fix to ieff to be single phase for qloss ca…

…lculation

src/PowerModelsGMD.jl

@@ -39,6 +39,9 @@ const _PMGMD = PowerModelsGMD
     import CSV
     using DataFrames
 
+    import Graphs
+    const _Gph = Graphs
+
     # Add core functions:
     include("core/base.jl")
     include("core/base_staged.jl")
@@ -79,6 +82,7 @@ const _PMGMD = PowerModelsGMD
     include("prob/gmd_mld.jl")
     include("prob/gmd_ots.jl")
     include("prob/gmd_ots_staged.jl")
+    include("prob/gmd_bounds.jl")
 
     include("prob/gmd_opf_staged.jl")
     include("prob/gmd_mld_staged.jl")

src/core/constraint.jl

@@ -195,7 +195,7 @@ end
 
 "CONSTRAINT: nodal power balance with gmd, shunts, and constant power factor load shedding"
 function constraint_power_balance_gmd_shunt_ls(pm::_PM.AbstractWConvexModels, n::Int, i::Int, bus_arcs, bus_arcs_dc, bus_arcs_sw, bus_gens, bus_storage, bus_pd, bus_qd, bus_gs, bus_bs)
-
+    println(ooo)
     w = _PM.var(pm, n, :w, i)
     p = get(_PM.var(pm, n), :p, Dict()); _PM._check_var_keys(p, bus_arcs, "active power", "branch")
     q = get(_PM.var(pm, n), :q, Dict()); _PM._check_var_keys(q, bus_arcs, "reactive power", "branch")
@@ -318,31 +318,6 @@ function constraint_qloss(pm::_PM.AbstractPowerModel, n::Int, k, i, j, baseMVA,
     end
 end
 
-"CONSTRAINT: qloss calculcated from dc current only"
-function constraint_qloss(pm::_PM.AbstractPowerModel, n::Int, k, i, j, baseMVA, branchMVA, busKV, K)
-    branch    = _PM.ref(pm, n, :branch, k)
-
-    qloss = _PM.var(pm, n, :qloss)
-    vm    = _PM.var(pm, n, :vm)[i]
-    ieff = _PM.var(pm, n, :i_dc_mag, k)
-
-
-
-    if branch["type"] == "xfmr"
-        JuMP.@constraint(pm.model,
-            qloss[(k,i,j)] == K * ieff * vm
-        )
-        # qloss[(k,i,j)] == 200.0
-    else
-        JuMP.@constraint(pm.model,
-            qloss[(k,i,j)] == 0.0
-        )
-    end
-        # Use this if we implement piecewise K
-        # (pm.data["baseMVA"]) / branchMVA ) * (K * vm * ieff) / (3.0 * branchMVA)
-        # (K * vm * ieff) / (3.0 * baseMVA)
-end
-
 
 "CONSTRAINT: qloss calculcated from ac voltage and constant ieff"
 function constraint_qloss_constant_ieff(pm::_PM.AbstractPowerModel, n::Int, k, i, j, baseMVA, K, ieff)
@@ -408,3 +383,11 @@ function constraint_dc_power_balance_ne_blocker(pm::_PM.AbstractPowerModel, n::I
     end
 
 end
+
+
+function constraint_gmd_connections(pm::_PM.AbstractPowerModel, n::Int, connections)
+    z = [_PM.var(pm, n, :z_blocker)[sub] for sub in connections]
+    con = JuMP.@constraint(pm.model,
+        sum([1 - _PM.var(pm, n, :z_blocker)[sub] for sub in connections]) >= 1
+    )
+end

src/core/constraint_template.jl

@@ -301,16 +301,15 @@ function constraint_qloss(pm::_PM.AbstractPowerModel, k; nw::Int=nw_id_default)
 
     branch    = _PM.ref(pm, nw, :branch, k)
     baseMVA   = _PM.ref(pm, :baseMVA)
-    branchMVA    = branch["baseMVA"]
+
     i         = branch["hi_bus"]
     j         = branch["lo_bus"]
 
     bus       = _PM.ref(pm, nw, :bus, i)
-    busKV     = bus["base_kv"]
 
     K         = calc_branch_K(pm,k;nw=nw)
 
-    constraint_qloss(pm, nw, k, i, j, baseMVA, branchMVA, busKV, K)
+    constraint_qloss(pm, nw, k, i, j, baseMVA, K)
 
 end
 
@@ -376,3 +375,36 @@ function constraint_dc_power_balance_ne_blocker(pm::_PM.AbstractPowerModel, i::I
     end
 
 end
+
+
+"CONSTRAINT: nodal current balance for dc circuits with GIC blockers and shunts"
+function constraint_dc_kcl_ne_blocker(pm::_PM.AbstractPowerModel, i::Int; nw::Int=nw_id_default)
+
+    dc_expr = pm.model.ext[:nw][nw][:dc_expr]
+    gmd_bus = _PM.ref(pm, nw, :gmd_bus, i)
+    gmd_bus_arcs = _PM.ref(pm, nw, :gmd_bus_arcs, i)
+    ne_blockers = get(_PM.ref(pm,nw,:gmd_bus_ne_blockers),i, Dict())
+    blockers = get(_PM.ref(pm,nw,:gmd_bus_blockers),i,Dict())
+
+    gs = gmd_bus["g_gnd"]
+    blocker_status = length(blockers) > 0 ? 1 : 0
+
+    if blocker_status == 0 && length(ne_blockers) > 0
+        if (length(ne_blockers) > 1)
+            Memento.warn(_LOGGER, "Bus ", i, " has more than one expansion blocker defined for it. Only using one of them")
+        end
+
+        constraint_dc_kcl_ne_blocker(pm, nw, i, ne_blockers[1], dc_expr, gmd_bus_arcs, gs)
+    else
+        constraint_dc_kcl(pm, nw, i, dc_expr, gmd_bus_arcs, gs, blocker_status)
+    end
+end
+
+
+"CONSTRAINT: that ensures that the connection does not become ungrounded"
+function constraint_gmd_connections(pm::_PM.AbstractPowerModel, i::Int; nw::Int=nw_id_default)
+    connections = _PM.ref(pm, nw, :gmd_connections, i)
+    if length(connections) >= 1 
+        constraint_gmd_connections(pm, nw, connections)
+    end
+end

src/core/data.jl

@@ -10,13 +10,34 @@
 
 "FUNCTION: calculate the minimum dc voltage at a gmd bus "
 function calc_min_dc_voltage(pm::_PM.AbstractPowerModel, i; nw::Int=pm.cnw)
-    return -1e6
+    gmd_bus = _PM.ref(pm, nw, :gmd_bus, i)
+    if haskey(gmd_bus, "vmin")
+        return gmd_bus["vmin"]
+    else
+        return -1e6
+    end
 end
 
 
 "FUNCTION: calculate the maximum dc voltage at a gmd bus "
 function calc_max_dc_voltage(pm::_PM.AbstractPowerModel, i; nw::Int=pm.cnw)
-    return 1e6
+    gmd_bus = _PM.ref(pm, nw, :gmd_bus, i)
+    if haskey(gmd_bus, "vmax")
+        return gmd_bus["vmax"]
+    else
+        return 1e6
+    end
+end
+
+
+"FUNCTION: calculate the max q loss of xfmr in branch "
+function calc_max_q_loss(pm::_PM.AbstractPowerModel, l; nw::Int=pm.cnw)
+    branch = _PM.ref(pm, nw, :branch, l)
+    if haskey(branch, "qloss_max")
+        return branch["qloss_max"] 
+    else
+        return 1e6
+    end
 end
 
 
@@ -100,12 +121,7 @@ function calc_ieff_current_mag_gwye_delta_xf(branch, case::Dict{String,Any}, sol
         if haskey(branch,"hi_3w_branch")
             return 0.0
         else
-            if haskey(solution["gmd_branch"]["$khi"], "dcf")
-                return abs(solution["gmd_branch"]["$khi"]["dcf"])
-            else
-                Memento.warn(_LOGGER, "Gwye-delta transformers $k doesn't have high-side dcf, skipping calculation of ieff")
-                return 0.0
-            end
+            return abs(solution["gmd_branch"]["$khi"]["dcf"])
         end
     end
 end
@@ -445,7 +461,7 @@ function calc_dc_current_line(branch, solution)
     g = 1 / branch["br_r"]
     vf = solution["gmd_bus"]["$nf"]["gmd_vdc"]
     vt = solution["gmd_bus"]["$nt"]["gmd_vdc"]
-    return g / 3 * (branch["br_v"] + (vf - vt))
+    return g * (branch["br_v"] + (vf - vt))
 end
 
 
@@ -455,7 +471,7 @@ function calc_dc_current_xfmr(branch, solution)
     g = 1 / branch["br_r"]
     vf = solution["gmd_bus"]["$nf"]["gmd_vdc"]
     vt = solution["gmd_bus"]["$nt"]["gmd_vdc"]
-    return g / 3 * (vf - vt)
+    return g * (vf - vt)
 end
 
 
@@ -472,14 +488,14 @@ end
 "FUNCTION: calculate the maximum DC current on a branch"
 function calc_dc_mag_max(pm::_PM.AbstractPowerModel, i; nw::Int=pm.cnw)
     branch = _PM.ref(pm, nw, :branch, i)
+    ibase = calc_branch_ibase(pm, i, nw=nw)
 
-    ac_max = -Inf
-    for l in _PM.ids(pm, nw, :branch)
-        ac_max = max(calc_ac_current_mag_max(pm, l, nw=nw), ac_max)
+    if haskey(branch, "qloss_max")
+        dc_mag_max = branch["qloss_max"]/(1.1*branch["gmd_k"])*ibase
+    else
+        dc_mag_max = 1e6
     end
 
-    ibase = calc_branch_ibase(pm, i, nw=nw)
-    dc_mag_max = 2 * ac_max * ibase
 
     if dc_mag_max < 0
         Memento.warn(_LOGGER, "DC current max for branch $i has been calculated as < 0. This will cause many things to break")
@@ -490,17 +506,17 @@ function calc_dc_mag_max(pm::_PM.AbstractPowerModel, i; nw::Int=pm.cnw)
 end
 
 
-"FUNCTION: calculate the ibase for a branch"
+"FUNCTION: calculate the ibase for a branch single phase"
 function calc_branch_ibase(pm::_PM.AbstractPowerModel, i; nw::Int=pm.cnw)
     branch = _PM.ref(pm, nw, :branch, i)
     bus = _PM.ref(pm, nw, :bus, branch["hi_bus"])
 
-    return branch["baseMVA"] * 1000.0 * sqrt(2.0) / (bus["base_kv"] * sqrt(3.0))
+    return branch["baseMVA"] * 1000.0 * sqrt(2.0) / (bus["base_kv"] * sqrt(3.0)) * 3
 end
 
 
 # ===   CALCULATIONS FOR QLOSS VARIABLES   === #
-
+"FUNCTION: returns qloss value in 3phase MVA"
 function calc_branch_K(pm::_PM.AbstractPowerModel, i; nw::Int=pm.cnw)
      branch = _PM.ref(pm, nw, :branch, i)
      ibase = calc_branch_ibase(pm,i;nw=nw)
@@ -509,8 +525,16 @@ function calc_branch_K(pm::_PM.AbstractPowerModel, i; nw::Int=pm.cnw)
 return haskey(branch, "gmd_k") ? (branch["gmd_k"] * branch["baseMVA"]) / (ibase) : 0.0
 end
 
+# ===   CALCULATIONS FOR QLOSS VARIABLES   === #
 
-"FUNCTION: calculate qloss
+function calc_branch_K_pu(pm::_PM.AbstractPowerModel, i; nw::Int=pm.cnw)
+    branch = _PM.ref(pm, nw, :branch, i)
+    ibase = calc_branch_ibase(pm,i;nw=nw)
+
+   return haskey(branch, "gmd_k") ? (branch["gmd_k"]) / (ibase) / 3 : 0.0
+end
+
+"FUNCTION: calculate qloss returns MVA 
     "
 function calc_qloss(branch::Dict{String,Any}, case::Dict{String,Any}, solution::Dict{String,Any})
     if branch["type"] == "xfmr"
@@ -556,7 +580,7 @@ function calc_qloss(branch::Dict{String,Any}, case::Dict{String,Any}, solution::
             i_dc_mag = abs(solution["ieff"]["$(branch["index"])"]) / ibase
 
             K = branch["gmd_k"]
-        return K * i_dc_mag * vm * case["baseMVA"]
+        return K * i_dc_mag * vm * case["baseMVA"] # needs to be checked based off branch baseMVA
         end
 
     end
@@ -804,7 +828,7 @@ function adjust_gmd_phasing!(result)
 
     gmd_branches = result["solution"]["gmd_branch"]
     for branch in values(gmd_branches)
-        branch["dcf"] /= 3
+        branch["dcf"] = branch["dcf"] / 3
     end
 
     return result
@@ -1017,3 +1041,62 @@ function add_gmd_3w_branch!(data::Dict{String,<:Any})
         end
     end
 end
+
+
+function get_connected_components(data::Dict{String,<:Any})
+    g = build_adjacency_matrix(data)
+    c = _Gph.SimpleGraph(g)
+    return _Gph.connected_components(c)
+end
+
+
+function filter_gmd_ne_blockers!(data::Dict{String,<:Any})
+    connections = []
+    for connection in data["connected_components"]
+        if length(connection) > 1
+            for conn in connection
+                append!(connections, conn)
+            end
+        end
+    end
+    gmd_ne_blocker = Dict{String, Any}()
+    for (i, gmd_blocker) in data["gmd_ne_blocker"]
+        if gmd_blocker["index"] in connections
+            gmd_ne_blocker[i] = gmd_blocker
+        end
+    end
+    data["gmd_ne_blocker"] = gmd_ne_blocker
+end
+
+
+function update_cost_multiplier!(data::Dict{String,<:Any})
+    subs = Dict{Int,Any}()
+    for (i, branch) in data["branch"]
+        if branch["type"] == "xfmr"
+            sub = []
+            branch_keys = ["gmd_br_hi", "gmd_br_lo", "gmd_br_series", "gmd_br_common"]
+            for branch_key in branch_keys
+                if branch[branch_key] != -1
+                    gmd_branch = data["gmd_branch"]["$(branch[branch_key])"]
+                    f_bus = gmd_branch["f_bus"]
+                    t_bus = gmd_branch["t_bus"]
+                    if data["gmd_bus"]["$f_bus"]["sub"] != -1 && !(data["gmd_bus"]["$f_bus"]["sub"] in sub)
+                        append!(sub, data["gmd_bus"]["$f_bus"]["sub"])
+                    end
+                end
+            end
+            for s in sub
+                if haskey(subs, s)
+                    subs[s] += 1
+                else
+                    subs[s] = 1
+                end
+            end
+        end
+    end
+    if haskey(data, "gmd_ne_blocker")
+        for (sub, m) in subs
+            data["gmd_ne_blocker"]["$sub"]["multipler"] = m
+        end
+    end
+end

src/core/objective.jl

@@ -15,3 +15,35 @@ function objective_blocker_placement_cost(pm::_PM.AbstractPowerModel)
     )
 
 end
+
+
+"OBJECTIVE: maximize the qloss of the ac model"
+function objective_max_qloss(pm::_PM.AbstractPowerModel, nw::Int=nw_id_default)
+    k = get(pm.setting,"qloss_branch",false)
+    branch = _PM.ref(pm, nw, :branch)[k]
+    i = branch["hi_bus"]
+    j = branch["lo_bus"]
+    return JuMP.@objective(pm.model, Max,
+        sum(_PM.var(pm, n, :qloss)[(k,i,j)]
+        for (n, nw_ref) in _PM.nws(pm))
+    )
+end
+
+
+"OBJECTIVE: max/min the dc voltage at the sub station"
+function objective_bound_gmd_bus_v(pm::_PM.AbstractPowerModel, nw::Int=nw_id_default)
+
+    bus = get(pm.setting,"gmd_bus",false)
+
+    if get(pm.setting,"max",false)
+        return JuMP.@objective(pm.model, Max,
+            sum(_PM.var(pm, n, :v_dc)[bus]
+            for (n, nw_ref) in _PM.nws(pm))
+        )
+    else
+        return JuMP.@objective(pm.model, Min,
+            sum(_PM.var(pm, n, :v_dc)[bus]
+            for (n, nw_ref) in _PM.nws(pm))
+        )
+    end
+end