Discrepancy BIM vs BFM (#286)

* found for discrepancy SOC BFM vs BIM
* Documentation fix
* Looser bounds on square current magnitude variable

Author: frederikgeth

docs/src/quickguide.md

@@ -86,15 +86,15 @@ This network data can be modified in the same way as the previous Matpower `.m`
 ## Inspecting AC and DC branch flow results
 The flow AC and DC branch results are not written to the result by default. To inspect the flow results, pass a settings Dict
 ```julia
-result = run_opf("case3_dc.m", ACPPowerModel, IpoptSolver(), setting = Dict("output" => Dict("line_flows" => true)))
+result = run_opf("case3_dc.m", ACPPowerModel, IpoptSolver(), setting = Dict("output" => Dict("branch_flows" => true)))
 result["solution"]["dcline"]["1"]
 result["solution"]["branch"]["2"]
 ```
 
-The losses of a AC or DC branch can be derived:
+The losses of an AC or DC branch can be derived:
 ```julia
-loss_ac =  Dict(name => data["p_to"]+data["p_from"] for (name, data) in result["solution"]["branch"])
-loss_dc =  Dict(name => data["p_to"]+data["p_from"] for (name, data) in result["solution"]["dcline"])
+loss_ac =  Dict(name => data["pt"]+data["pf"] for (name, data) in result["solution"]["branch"])
+loss_dc =  Dict(name => data["pt"]+data["pf"] for (name, data) in result["solution"]["dcline"])
 ```
 
 

src/core/data.jl

@@ -44,30 +44,57 @@ function calc_series_current_magnitude_bound(branches, buses)
         g_sh_to = branch["g_to"]
         b_sh_fr = branch["b_fr"]
         b_sh_to = branch["b_to"]
-        zmag_fr = abs(g_sh_fr + im*b_sh_fr)
-        zmag_to = abs(g_sh_to + im*b_sh_to)
+        r_s = branch["br_r"]
+        x_s = branch["br_x"]
 
         vmax_fr = bus_fr["vmax"]
         vmax_to = bus_fr["vmax"]
-        vmin_fr = bus_fr["vmin"]
-        vmin_to = bus_fr["vmin"]
 
         tap_fr = branch["tap"]
         tap_to = 1 # no transformer on to side, keeps expressions symmetric.
         smax = branch["rate_a"]
 
-        cmax_tot_fr = smax*tap_fr/vmin_fr
-        cmax_tot_to = smax*tap_to/vmin_to
+        cmax_p = (2*smax + abs(g_sh_fr)*vmax_fr^2 + abs(g_sh_to)*vmax_to^2)/(abs(r_s))
+        cmax_q = (2*smax + abs(b_sh_fr)*vmax_fr^2 + abs(b_sh_to)*vmax_to^2)/(abs(x_s))
 
-        cmax_sh_fr = zmag_fr * vmax_fr
-        cmax_sh_to = zmag_to * vmax_to
-
-        cmax[key] = max(cmax_tot_fr + cmax_sh_fr, cmax_tot_to + cmax_sh_to)
+        cmax[key] = min(cmax_p, cmax_q)
     end
     return cmax
 end
 
 
+""
+function calc_series_active_power_bound(branches, buses)
+    pmax = Dict([(key, 0.0) for key in keys(branches)])
+    for (key, branch) in branches
+        bus_fr = buses[branch["f_bus"]]
+        g_sh_fr = branch["g_fr"]
+        vmax_fr = bus_fr["vmax"]
+        tap_fr = branch["tap"]
+        smax = branch["rate_a"]
+
+        pmax[key] = smax + abs(g_sh_fr) * (vmax_fr/tap_fr)^2
+    end
+    return pmax
+end
+
+
+""
+function calc_series_reactive_power_bound(branches, buses)
+    qmax = Dict([(key, 0.0) for key in keys(branches)])
+    for (key, branch) in branches
+        bus_fr = buses[branch["f_bus"]]
+        b_sh_fr = branch["g_fr"]
+        vmax_fr = bus_fr["vmax"]
+        tap_fr = branch["tap"]
+        smax = branch["rate_a"]
+
+        qmax[key] = smax + abs(b_sh_fr) * (vmax_fr/tap_fr)^2
+    end
+    return qmax
+end
+
+
 ""
 function calc_theta_delta_bounds(data::Dict{String,Any})
     bus_count = length(data["bus"])

src/form/df.jl

@@ -25,11 +25,14 @@ end
 
 
 function variable_active_branch_series_flow(pm::GenericPowerModel, n::Int=pm.cnw; bounded = true)
+    branches = pm.ref[:nw][n][:branch]
+    buses = pm.ref[:nw][n][:bus]
+    pmax = calc_series_active_power_bound(branches, buses)
     if bounded
         pm.var[:nw][n][:p_s] = @variable(pm.model,
             [(l,i,j) in pm.ref[:nw][n][:arcs_from]], basename="$(n)_p_s",
-            lowerbound = -pm.ref[:nw][n][:branch][l]["rate_a"],
-            upperbound =  pm.ref[:nw][n][:branch][l]["rate_a"],
+            lowerbound = -pmax[l],
+            upperbound =  pmax[l],
             start = getstart(pm.ref[:nw][n][:branch], l, "p_start")
         )
     else
@@ -42,11 +45,14 @@ end
 
 
 function variable_reactive_branch_series_flow(pm::GenericPowerModel, n::Int=pm.cnw; bounded = true)
+    branches = pm.ref[:nw][n][:branch]
+    buses = pm.ref[:nw][n][:bus]
+    qmax = calc_series_reactive_power_bound(branches, buses)
     if bounded
         pm.var[:nw][n][:q_s] = @variable(pm.model,
             [(l,i,j) in pm.ref[:nw][n][:arcs_from]], basename="$(n)_q_s",
-            lowerbound = -pm.ref[:nw][n][:branch][l]["rate_a"],
-            upperbound =  pm.ref[:nw][n][:branch][l]["rate_a"],
+            lowerbound = -qmax[l],
+            upperbound =  qmax[l],
             start = getstart(pm.ref[:nw][n][:branch], l, "q_start")
         )
     else
@@ -96,11 +102,11 @@ function constraint_voltage_magnitude_difference(pm::GenericPowerModel{T}, n::In
     ccm =  pm.var[:nw][n][:ccm][i]
 
     #define series flow expressions to simplify Ohm's law
-    p_fr_s = p_fr - g_sh_fr*(w_fr/tm^2)
-    q_fr_s = q_fr + b_sh_fr*(w_fr/tm^2)
+    p_s_fr = p_fr - g_sh_fr*(w_fr/tm^2)
+    q_s_fr = q_fr + b_sh_fr*(w_fr/tm^2)
 
     #KVL over the line:
-    @constraint(pm.model, w_to == (w_fr/tm^2) - 2*(r*p_fr_s + x*q_fr_s) + (r^2 + x^2)*ccm)
+    @constraint(pm.model, w_to == (w_fr/tm^2) - 2*(r*p_s_fr + x*q_s_fr) + (r^2 + x^2)*ccm)
 
 end
 

test/opf.jl

@@ -294,7 +294,7 @@ end
         result = run_opf_bf("../test/data/matpower/case5_gap.m", SOCDFPowerModel, ipopt_solver)
 
         @test result["status"] == :LocalOptimal
-        @test isapprox(result["objective"], -27660.4; atol = 1e0)
+        @test isapprox(result["objective"], -28238.0; atol = 1e0)
     end
     @testset "5-bus with pwl costs" begin
         result = run_opf_bf("../test/data/matpower/case5_pwlc.m", SOCDFPowerModel, ipopt_solver)