Debugging and validation: work in progress

[frederikgeth]

Warning: long post

I have some observations, after playing with the current implementation in a fork. You can run this script if you do a fix for https://github.com/lanl-ansi/ThreePhasePowerModels.jl/issues/45 in the latest master of PowerModels.jl, e.g.

function calc_branch_y(branch::Dict{String,Any})
    y = (branch["br_r"].values + im*branch["br_x"].values)^-1
    g = real(y)
    b = imag(y)
    return g, b
end

I largely followed [1], to develop the SDP BFM and the 'simplified' LP BFM. I also added a SOC formulation based on the SOC relaxation of the PSD constraints. I don't trust the SDP and SOC (both as NLP and as conic, for debugging) formulations at the moment. Nevertheless, the NLP solved in the script is consistent (losses, flows, voltages) with OpenDSS to <1e-5 and the LP correlates strongly, so I sort of trust those. The equations are developed using JuMP's vectorised syntax, therefore I also added a new problem type tp_opf_bf_mat.

Before I attempt at converting this to a PR, it would be nice to exchange ideas on the following:

• The Hermitian matrix variables for current ccm_re + j ccm_im and voltage w_re + j w_im re-use the variables PMs.variable_voltage_magnitude_sqr,PMs.variable_current_magnitude_sqr on the diagonal.
• New vector variables ccmr, wr and ccmi, wi are defined to represent the off-diagonal real and imaginary elements of ccm_re + j ccm_im and w_re + j w_im respectively.
• The full matrix variable for variable_branch_flow re-uses PMs.variable_active_branch_flow for the diagonal and new variables for upper and lower triangle elements (with suffix _ut, _lt).
• Functions in multiphase.jl allow for constructing matrices from these vectors, and extracting the lower/upper triangles from such matrices.
• The base kV read in from OpenDSS refers to the line-line voltage (e.g. 400 V) whereas the voltage variables actually are line-neutral, so I had to adapt the basekv
• Somehow the impedance was off as well (don't have this one figured out).
• Where to put the algorithm ('Algorithm 2' in [1], crappily implemented in post_process_voltage_current in the script) to derive the natural variables

Todos:

• [x] separate model variants of ref_bus: e.g. bus with fixed voltage magnitude, bus with fixed angle, or both. OpenDSS fixes the voltage at the sourcebus based on the pu setting of the circuit code.
• [x] implement tighter bounds for off-diagonal elements of the matrix variables
• [x] distinguish constraints defined over phases from defined over conductors (going towards 4x4 and 5x5)
• [x] work out constraints in a scalar manner / ensure compatibility with tp_opf_bf
• [x] remove debug code
• [x] add unit tests
• [x] understand infeasibility of SDP/SOC
• [x] documentation
• [x] remove .values after fix for https://github.com/lanl-ansi/PowerModels.jl/issues/313
• [x] remove constraint_tp_valid_inquality_branch and constraint_tp_valid_inquality_bus
• [ ] merge all bf formulations in one file
• [x] trace(p_fr) + trace(p_to) >= 0 vs diag(p_fr) + diag(p_to) >= 0

Comments welcome!

-- [1] Gan, L., & Low, S. H. (2014). Convex relaxations and linear approximation for optimal power flow in multiphase radial networks. PSSC

Reference results OpenDSS for case3_unbalanced:

Bus   (node ref)  Node       V (kV)    Angle    p.u.   Base kV
SOURCEBUS .. 1      0.22999 /_    0.0    0.9959     0.400        1-2    0.39836 /_   30.0    0.9959
   -         2      0.22999 /_ -120.0    0.9959     0.400        2-3    0.39836 /_  -90.0    0.9959
   -         3      0.22999 /_  120.0    0.9959     0.400        3-1    0.39836 /_  150.0    0.9959
PRIMARY .... 1      0.22654 /_   -0.2   0.98094     0.400        1-2    0.39384 /_   30.0    0.9846
   -         2      0.22848 /_ -120.1   0.98936     0.400        2-3    0.39481 /_  -90.0   0.98702
   -         3      0.22795 /_  120.1   0.98704     0.400        3-1    0.39428 /_  149.8   0.98571
LOADBUS .... 1      0.22252 /_   -0.5   0.96355     0.400        1-2    0.38859 /_   29.9   0.97148
   -         2      0.22673 /_ -120.2   0.98176     0.400        2-3    0.39068 /_  -90.1   0.97671
   -         3      0.22558 /_  120.3   0.97678     0.400        3-1    0.38954 /_  149.7   0.97385

   Bus       Phase     kW     +j   kvar           kVA           PF
ELEMENT = "Line.OHLINE"   
SOURCEBUS       1     9.2757 +j     3.1776        9.8048      0.9460
SOURCEBUS       2     6.0735 +j     3.0673        6.8041      0.8926
SOURCEBUS       3      6.132 +j     3.0278        6.8388      0.8967
 TERMINAL TOTAL       21.481 +j     9.2726        23.397      0.9181
ELEMENT = "Line.OHLINE"   
PRIMARY         1    -9.1485 +j    -3.0949        9.6578      0.9473
PRIMARY         2    -6.0396 +j     -3.036        6.7598      0.8935
PRIMARY         3     -6.071 +j    -3.0147        6.7783      0.8957
 TERMINAL TOTAL      -21.259 +j    -9.1456        23.143      0.9186

ELEMENT = "Line.QUAD"     
PRIMARY         1     9.1485 +j     3.0949        9.6578      0.9473
PRIMARY         2     6.0396 +j      3.036        6.7598      0.8935
PRIMARY         3      6.071 +j     3.0147        6.7783      0.8957
 TERMINAL TOTAL       21.259 +j     9.1456        23.143      0.9186
ELEMENT = "Line.QUAD"     
LOADBUS         1         -9 +j         -3        9.4868      0.9487
LOADBUS         2         -6 +j         -3        6.7082      0.8944
LOADBUS         3         -6 +j         -3        6.7082      0.8944
 TERMINAL TOTAL          -21 +j         -9        22.847      0.9191

[ccoffrin]

Sounds like some great progress. I suggest that we setup a call to discuss any details early next week.

Some initial questions

• The current code implements an SOC relaxation under SOCWRPowerModel, is your goal to have another SOC relaxation based on the BFM formulation? Can you use the current SOC as a refrence for your current-based version?
• Is there a notable difference between LP BFM and LinDistFlow?

[frederikgeth]

• I implemented Gan and Low's unbalanced BFM formulation referenced above, as they note that the SDP relaxation of the unbalanced BIM formulation does not have very nice numerical properties. They don't define a SOC relaxation of unbalanced BFM though.
• To debug the SDP formulation, I added some constraints to define valid inequalities (e.g. explicit PSDness of the current and voltage matrix variables) and performed the SDP->SOC convex relaxation (all 2x2 matrices formed by taking any two diagonal elements ii and jj, and then also ij and ji should also be PSD, what can be cast as a rotated SOC, per [2]). I added support for solving this as an NLP and as a true conic problem, because I ran into numerical issues.
• The LP BFM formulation of Gan and Low isn't identical to the LinDistFlow formulation, as LinDistFlow only defines equality on the diagonal in this constraint. But I'm not sure they actually are different (due to the properties of the matrices involved), this is to be tested.

[2] Kim, Kojima and Yamashita, Second Order Cone Programming Relaxation of a Positive Semidefinite Constraint, Optimization Methods and Software, 2003

[frederikgeth]

Yeah, and a call next week is a great idea!

[frederikgeth]

Some helper functions to consider for the line impedance matrices here.

The license of the test data is to be discussed before making it a part of this package, see:

• https://dspace.lboro.ac.uk/dspace-jspui/handle/2134/15544
• https://creativecommons.org/licenses/by-nc/3.0/

[frederikgeth]

I will close this now, major issues are addressed or taken up in more targeted open issues.