Extending OBBT in PowerModels

I want to extend the current version of OBBT approach in PowerModels by tightening “Voltage Magnitude Differences” besides the “Voltage Magnitude” and “Voltage Angle Dffference” variables. In this connection, I have modified branch matrix in Matpower test cases by adding six new columns (22-27). Note that I have reserved columns (14-21) for post_pf or post_opf results (i.e., pf, qf, etc.). I have also modified “Matpowe.jl” script to take care of the new six added columns to the branch data matrix in Matpower.
Regarding the OBBT script itself, I have added some codes for tightening “Voltage Magnitude Difference”. It was not complicated because tightening “Voltage Magnitude Differences” has a similar procedure as tightening “Voltage Angle Difference”. However, I have some difficulties with updating data after each iteration in OBBT. I do know there is a script in InfrastructureModels (i.e., data.jl) which is used to update data, but I got confused when I wanted to modify that.
I want to write my code in the most general way in order to contribute to PowerModels later. Thus, I would really appreciate it if you help me to get this done. Any hint would be appreciated!

If your goal is to do OBBT, you will need to extend one of the PowerModels formulations to include an explicit variable for the voltage magnitude difference and a linking constraint that connects that quantity back to the original model (otherwise it brings no additional value to the OBBT algorithm).

I think the simplest route is to make a new PowerModel type that is a sub-type of QCWRPowerModel or QCWRTriPowerModel . With this new type you will want to copy and extend the functions for variable_voltage (add the new voltage magnitude difference variables), constraint_voltage (add the new linking constraint).

Regarding the extension of the Matpower data parser. You should not need to hack any of the current codebase to do this. PowerModels supports ad-hoc extensions of the Matpower format, see this part of the docs. In short, add a new matrix to the Matpower file with the correct naming convention (i.e. mpc.branch_) and the same number of rows as the branch matrix and PowerModels’ Matpower parser will add the new data into the branch dictionary objects for you.

That said, you may prefer to simply add the voltage magnitude difference bounds to the PowerModels data as a post-processing step, after parsing the Matpower file. This would allow your method to work with both Matpower and PTI data files.

Thank you so much for the response! I modified Matpower test cases by adding a new matrix in order to take care of six new parameters that I wanted to be added to mpc.branch. I have already constructed a new type by adding the corresponding variables and constraints related to “Voltage magnitude difference” variables. The new type is a subtype of the QCWRTriPowerModel.

I have modified base.jl by adding “vdiff_min” to buspairs as parameters. It was straight forward since I have followed the same steps that had been taken for “angmin” and “angmax” parameters.

I also had to modified data.jl in order to take care of the sign of added parameters for flipped branches. It was also straight forward because it imitated what had been done for “angmax” and “angmin”.

I also added new variables corresponds to “voltage magnitude difference” to the new type. Following are the variables that I have defined.


function variable_voltage_magnitude_difference_neghboring_buses(pm::GenericPowerModel{QCWRTriForm}; nw::Int=pm.cnw, cnd::Int=pm.ccnd)

##buspairs = ref(pm, nw, :buspairs)

var(pm, nw, cnd)[:vd] = @variable(pm.model,

[bp in ids(pm, nw, :buspairs)], basename="$(nw)_$(cnd)_vd",

##[bp in keys(buspairs)], basename="$(nw)_$(cnd)_vd",

lowerbound = ref(pm, nw, :buspairs, bp, "vdiff_min", cnd),

upperbound = ref(pm, nw, :buspairs, bp, "vdiff_max", cnd),

start = getval(ref(pm, nw, :buspairs, bp), "vd_start", cnd)

)

end

function variable_voltage_difference_product_voltage_from(pm::GenericPowerModel{QCWRTriForm}; nw::Int=pm.cnw, cnd::Int=pm.ccnd)

buspairs = ref(pm, nw, :buspairs)

var(pm, nw, cnd)[:vvv_fr] = @variable(pm.model,

# [bp in ids(pm, nw, :buspairs)], basename="$(nw)_$(cnd)_vvv_fr",

[bp in keys(buspairs)], basename="$(nw)_$(cnd)_vvv_fr",

lowerbound = buspairs[bp]["vm_fr_min"][cnd]*buspairs[bp]["vdiff_min"][cnd],

upperbound = buspairs[bp]["vm_fr_max"][cnd]*buspairs[bp]["vdiff_max"][cnd],

start = getval(ref(pm, nw, :buspairs, bp), "vvv_fr_start", cnd)

)

 end

function variable_voltage_difference_product_voltage_to(pm::GenericPowerModel{QCWRTriForm}; nw::Int=pm.cnw, cnd::Int=pm.ccnd)

buspairs = ref(pm, nw, :buspairs)

var(pm, nw, cnd)[:vvv_to] = @variable(pm.model,

## [bp in ids(pm, nw, :buspairs)], basename="$(nw)_$(cnd)_vvv_to",

[bp in keys(buspairs)], basename="$(nw)_$(cnd)_vvv_to",

lowerbound = buspairs[bp]["vm_to_min"][cnd]*buspairs[bp]["vdiff_min"][cnd],

upperbound = buspairs[bp]["vm_to_max"][cnd]*buspairs[bp]["vdiff_max"][cnd],

start = getval(ref(pm, nw, :buspairs, bp), "vvv_to_start", cnd)

)

end

function variable_voltage_squared_product(pm::GenericPowerModel{QCWRTriForm}; nw::Int=pm.cnw, cnd::Int=pm.ccnd)

buspairs = ref(pm, nw, :buspairs)

var(pm, nw, cnd)[:w_squared] = @variable(pm.model,

#[bp in ids(pm, nw, :buspairs)], basename="$(nw)_$(cnd)_w_squared",

[bp in keys(buspairs)], basename="$(nw)_$(cnd)_w_squared",

lowerbound = buspairs[bp]["vm_fr_min"][cnd]^2*buspairs[bp]["vm_to_min"][cnd]^2,

upperbound = buspairs[bp]["vm_fr_max"][cnd]^2*buspairs[bp]["vm_to_max"][cnd]^2,

start = getval(ref(pm, nw, :buspairs, bp), "w_squared", cnd)

)

 end

function variable_voltage_magnitude_diff_sqr(pm::GenericPowerModel{QCWRTriForm}; nw::Int=pm.cnw, cnd::Int=pm.ccnd)

buspairs = ref(pm, nw, :buspairs)

var(pm, nw, cnd)[:vd_squared] = @variable(pm.model,

## [bp in ids(pm, nw, :buspairs)], basename="$(nw)_$(cnd)_vd_squared",

[bp in keys(buspairs)], basename="$(nw)_$(cnd)_vd_squared",

lowerbound = min(buspairs[bp]["vdiff_min"][cnd]^2,buspairs[bp]["vdiff_max"][cnd]^2),

upperbound = max(buspairs[bp]["vdiff_min"][cnd]^2,buspairs[bp]["vdiff_max"][cnd]^2),

start = getval(ref(pm, nw, :buspairs, bp), "vd_squared", cnd)

)

end

As you mentioned I added these variables to the “variable_voltage”. I have also defined few constraints related to voltage magnitude difference variables as follows


function constraint_voltage_magnitude_difference_neghboring_buses(pm::GenericPowerModel{QCWRTriForm}, n::Int, c::Int, f_idx, vdiff_min, vdiff_max)

i, f_bus, t_bus = f_idx

vm_fr = var(pm, n, c, :vm, f_bus)

vm_to = var(pm, n, c, :vm, t_bus)

@constraint(pm.model, vm_fr - vm_to <= vdiff_max)

@constraint(pm.model, vm_fr - vm_to >= vdiff_min)

end

function constraint_voltage_magnitude_diff_sqr_equality(pm::GenericPowerModel{QCWRTriForm}, n::Int, c::Int, f_bus, t_bus, arc_from, tm)

w_fr = var(pm, n, c, :w, f_bus)

w_to = var(pm, n, c, :w, t_bus)

vv = var(pm, n, c, :vv, (f_bus, t_bus))

vd = var(pm, n, c, :vd, (f_bus, t_bus))

@constraint(pm.model, vv == (w_fr + w_to - vd_squared)/2)

end

function constraint_voltage_magnitude_diff_sqr_inequality(pm::GenericPowerModel{QCWRTriForm}, n::Int, c::Int, f_bus, t_bus)

w_fr = var(pm, n, c, :w, f_bus)

w_to = var(pm, n, c, :w, t_bus)

vv = var(pm, n, c, :vv, (f_bus, t_bus))

vd = var(pm, n, c, :vd, (f_bus, t_bus))

@constraint(pm.model, vd_squared <= w_fr - 2*vv + w_to)

end

function constraint_voltage_magnitude_diff_link(pm::GenericPowerModel{QCWRTriForm}, n::Int, c::Int, f_bus, t_bus, g, b, b_fr)

p_fr = var(pm, n, c, :p, f_idx)

q_fr = var(pm, n, c, :q, f_idx)

q_to = var(pm, n, c, :q, t_idx)

p_to = var(pm, n, c, :p, t_idx)

vd = var(pm, n, c, :vd, (f_bus, t_bus))

vvv_fr = var(pm, n, c, :vvv_fr, f_bus)

 vvv_to = var(pm, n, c, :vvv_to, t_bus)

 @constraint(pm.model, vvv_fr + vvv_to == (1/(g^2+b^2+b*b_fr))*(g*(p_fr-p_to)-b*(q_fr-q_to)))

end

# Following constraints are Dmitry's constraints in the ShareLatex file

 function constraint_voltage_magnitude_diff_dim1(pm::GenericPowerModel{QCWRTriForm}, n::Int, c::Int, f_idx,alpha_fr, alpha_to, beta_fr, beta_to)

i, f_bus, t_bus = f_idx

vm_fr = var(pm, n, c, :vm, f_bus)

vm_to = var(pm, n, c, :vm, t_bus)

@constraint(pm.model, vm_to - alpha_fr*vm_fr <= beta_fr)

@constraint(pm.model, vm_to - alpha_fr*vm_fr >= -beta_fr)

end

function constraint_voltage_magnitude_diff_dim2(pm::GenericPowerModel{QCWRTriForm}, n::Int, c::Int, f_idx, alpha_fr, alpha_to, beta_fr, beta_to)

i, f_bus, t_bus = f_idx

vm_fr = var(pm, n, c, :vm, f_bus)

vm_to = var(pm, n, c, :vm, t_bus)

w_fr = var(pm, n, c, :w, f_bus)

w_to = var(pm, n, c, :w, t_bus)

vv = var(pm, n, c, :vv, (f_bus, t_bus))

@constraint(pm.model, w_to - 2*alpha_fr*vv + alpha_fr^2 *w_fr <= beta_fr^2)

end

function constraint_voltage_magnitude_diff_dim3(pm::GenericPowerModel{QCWRTriForm}, n::Int, c::Int, f_idx, alpha_fr, alpha_to, beta_fr, beta_to)

i, f_bus, t_bus = f_idx

vm_fr = var(pm, n, c, :vm, f_bus)

vm_to = var(pm, n, c, :vm, t_bus)

vv = var(pm, n, c, :vv, (f_bus, t_bus))

@constraint(pm.model, w_to + 2*alpha_fr*vv + alpha_fr^2 * w_fr <= 4*alpha_fr*vv + beta_fr^2)

End

When I put these constraint as function in “constraint_voltage” the code worked without any error, but the result of bound tightening was not correct. When I defined these constraints outside of the “constraint_voltage” I had some difficulties to call them in “constraint_voltage” because parameters such as “alpha_fr” were not recognized. I will really appreciate it if you help me to figure out the issue. Thanks in advance!

If you would like assistance at this level of detail, I think a PR to PowerModelsAnnex would be the best venue. Note in that case, your implementation will need to work without modifying the core of PowerModels. If you would like a feature added to PowerModels, post an issue to PowerModels.jl and we can discuss it there.

Thanks for the response! Sure, I will do that in the future. Just two quick questions, as you suggested I have modified the “test case.m” by adding an extra matrix in order to add my desired values. I have also extended the OBBT.jl to tighten the voltage difference as well. My first question is, does update_data function recognize the new variable (i.e., voltage difference bounds) or not? I just want to be sure the issue does not stem from that function.
My second question is about create_modifications function in OBBT.jl. There are two dictionaries have been defined (modifications[“bus”] = Dict{String,Any}()
and modifications[“branch”] = Dict{String,Any}()) I need the third dictionary to take care of modification to voltage magnitude differences. The problem is the voltage magnitude difference is also a dictionary of branches. Should I define the dictionary for modification of “voltage magnitude difference variable” based on the matrix that I have added to the test_case.m file?
I will appreciate your response! I have not touched the core of PowerModels in my new implementation.

does update_data function recognize the new variable (i.e., voltage difference bounds) or not?

Yes, update_data is designed to work with arbitrary dictionaries.

My second question is about create_modifications function in OBBT.jl …

You would want to add your list of new bounds to the function arguments and then modify the following line,

modifications["branch"][index] = Dict{String,Any}( "angmin" => td_lb[bp], "angmax" => td_ub[bp] )

To something along the lines of,

modifications["branch"][index] = Dict{String,Any}( "angmin" => td_lb[bp], "angmax" => td_ub[bp], "vmdmin" => vmd_lb[bp], "vmdmax" => vmd_ub[bp])

Thank you so much for the instructive response!