Reproducing power flow results

[antonhinneck]

To show the validity of my experiments, I want to use PowerModels.jl to provide reference solutions.
For my current work, I am specifically interested in the SOC approximation - corresponding to SOCWRPowerModel.

I had little success, however, consistently reproducing results.

Thus, I went on to trying to reproduce SOCWRPowerModel as exact as possible.
Here is the current model:
I used the notation in Jabr 2006, so wr = Rij, wi = Iij, w = ui.
For "pglib_opf_case3_lmbd.m" I included a comparison of the printed JuMP models of PowerModels.jl (left) and the model below (right).
In the model below pt[e] corresponds to p(e, fbus[e], tbus[e]), while pf[e] corresponds to p(e, tbus[e], fbus[e]).

While PowerModels.jl produces the following accurate result:
Dict{String, Any} with 3 entries:
"1" => Dict{String, Any}("qg"=>0.224946, "pg"=>0.950497)
"2" => Dict{String, Any}("qg"=>0.18509, "pg"=>2.24491)
"3" => Dict{String, Any}("qg"=>0.096942, "pg"=>0.0)

The model below produces this result:
Dict{String, Dict{String, Float64}} with 3 entries:
"1" => Dict("qg"=>-0.997569, "pg"=>0.756516)
"2" => Dict("qg"=>-0.514129, "pg"=>2.24656)
"3" => Dict("qg"=>-0.483644, "pg"=>1.0e-8)

I most recently

• added the start values,
• additional constraints to tighten the problem (lifted nonlinear cuts) and
• a lower tolerance of the solver,
  with no success.

Are there any glaring issues that I am just missing, or anything else PowerModels does under the hood?

If you have any suggestions, I would appreciate it a lot.

Best,

Anton

function soc_j07s(d; optimizer = optimizer_with_attributes(Ipopt.Optimizer, "tol" => 1e-6), obj = :inj)
    
    m = Model(optimizer)

    @variable(m, Iij[e = d.Branches], start = 0.0)
    @variable(m, Rij[e = d.Branches] >= 0.0, start = 0.0)
    @variable(m, ui[d.Buses] >= 0.0, start = 1.0)
    @variable(m, pf[e = d.Branches])
    @variable(m, pt[e = d.Branches])
    @variable(m, qf[e = d.Branches])
    @variable(m, qt[e = d.Branches])
    @variable(m, pg[g = d.Generators] >= 0.0)
    @variable(m, qg[d.Generators])

    @objective(m, Min, sum(d.c1[g] * pg[g] for g in d.Generators))

    if obj == :inj
        @objective(m, Min, sum(d.c1[g] * pg[g] for g in d.Generators))
    elseif obj == :j06
        @objective(m, Max, sum(Rij[e] for e in d.Branches))
    else
        error("Objective function: not defined ...")
    end

    @constraint(m, set_slack, ui[d.sb] == 1.0)

    @constraint(m, p_power_balance[n = d.Buses], sum(  pt[e] for e in d.Etn[n]) + sum(pf[e] for e in d.Efn[n]) + d.Pdn[n] == sum(pg[g] for g in d.Gn[n]) + d.Gsn[n])
    @constraint(m, q_power_balance[n = d.Buses], sum(  qt[e] for e in d.Etn[n]) + sum(qf[e] for e in d.Efn[n]) + d.Qdn[n] == sum(qg[g] for g in d.Gn[n]) + d.Bsn[n])

    @constraint(m, pf_def[e = d.Branches], pf[e] ==  ui[d.tbus[e]] * (d.g[e]+d.g_fr[e]) / d.tm[e]^2 
                                                    + (-d.g[e] * d.tr[e] + d.b[e] * d.ti[e]) * Rij[e] / d.tm[e]^2 
                                                    + (-d.b[e] * d.tr[e] - d.g[e] * d.ti[e]) * Iij[e] / d.tm[e]^2)
    @constraint(m, qf_def[e = d.Branches], qf[e] == -ui[d.tbus[e]] * (d.b[e]+d.b_fr[e]) / d.tm[e]^2 
        	                                        - (-d.b[e] * d.tr[e] - d.g[e] * d.ti[e]) * Rij[e] / d.tm[e]^2 
                                                    + (-d.g[e] * d.tr[e] + d.b[e] * d.ti[e]) * Iij[e] / d.tm[e]^2)

    @constraint(m, pt_def[e = d.Branches], pt[e] ==  ui[d.fbus[e]] * (d.g[e]+d.g_to[e])
                                                    + (-d.g[e] * d.tr[e] - d.b[e] * d.ti[e]) * Rij[e] / d.tm[e]^2 
                                                    + (-d.b[e] * d.tr[e] + d.g[e] * d.ti[e]) *-Iij[e] / d.tm[e]^2)
    @constraint(m, qt_def[e = d.Branches], qt[e] == -ui[d.fbus[e]] * (d.b[e]+d.b_to[e]) 
        	                                        - (-d.b[e] * d.tr[e] + d.g[e] * d.ti[e]) * Rij[e] / d.tm[e]^2 
                                                    + (-d.g[e] * d.tr[e] - d.b[e] * d.ti[e]) *-Iij[e] / d.tm[e]^2)

    @constraint(m, angle1[e = d.Branches], Iij[e] <= tan(d.amax[e])*Rij[e])
    @constraint(m, angle2[e = d.Branches], Iij[e] >= tan(d.amin[e])*Rij[e])

    if optimizer == Ipopt.Optimizer || optimizer == SCIP.Optimizer
        @NLconstraint(m, surf[e = d.Branches], ui[d.tbus[e]] * ui[d.fbus[e]] >= Rij[e]^2 + Iij[e]^2)
        @NLconstraint(m, lim_p_f[e = d.Branches], d.ra[e] * d.ra[e] >= pf[e]^2 + qf[e]^2)
        @NLconstraint(m, lim_p_t[e = d.Branches], d.ra[e] * d.ra[e] >= pt[e]^2 + qt[e]^2)
    end

    # Lifted Nonlinear Cuts
    #######################
    @expression(m, vf_ub[e=d.Branches], d.vmax[d.fbus[e]])
    @expression(m, vf_lb[e=d.Branches], d.vmin[d.fbus[e]])
    @expression(m, vt_ub[e=d.Branches], d.vmax[d.tbus[e]])
    @expression(m, vt_lb[e=d.Branches], d.vmin[d.tbus[e]])
    @expression(m, sf[e=d.Branches], d.vmin[d.fbus[e]] + d.vmax[d.fbus[e]])
    @expression(m, st[e=d.Branches], d.vmin[d.tbus[e]] + d.vmax[d.tbus[e]])
    @expression(m, phi1[e=d.Branches], (d.amax[e] + d.amin[e])/2)
    @expression(m, phi2[e=d.Branches], (d.amax[e] - d.amin[e])/2)
    @expression(m, wt[e=d.Branches], ui[d.tbus[e]])
    @expression(m, wf[e=d.Branches], ui[d.fbus[e]])

    @constraint(m, lifted_nonlinear_cuts_1[e = d.Branches], sf[e]*st[e]*(cos(phi1[e])*Rij[e]+sin(phi1[e])*Iij[e]) - vt_ub[e]*cos(phi2[e])*st[e]*wf[e] - vf_ub[e]*cos(phi2[e])*sf[e]*wt[e] >=  vf_ub[e]*vt_ub[e]*cos(phi2[e])*(vf_lb[e]*vt_lb[e] - vf_ub[e]*vt_ub[e]))
    @constraint(m, lifted_nonlinear_cuts_2[e = d.Branches], sf[e]*st[e]*(cos(phi1[e])*Rij[e]+sin(phi1[e])*Iij[e]) - vt_lb[e]*cos(phi2[e])*st[e]*wf[e] - vf_lb[e]*cos(phi2[e])*sf[e]*wt[e] >= -vf_lb[e]*vt_lb[e]*cos(phi2[e])*(vf_lb[e]*vt_lb[e] - vf_ub[e]*vt_ub[e]))
    
    ## WR, WI lower bounds
    ######################
    @constraint(m, wr_ub[e = d.Branches], Rij[e] <= d.vmax[d.fbus[e]]*d.vmax[d.tbus[e]]*cos(0))
    @constraint(m, wr_lb[e = d.Branches], Rij[e] >= d.vmin[d.fbus[e]]*d.vmin[d.tbus[e]]*cos(d.amax[e]))
    @constraint(m, wi_ub[e = d.Branches], Iij[e] <= d.vmax[d.fbus[e]]*d.vmax[d.tbus[e]]*sin(d.amax[e]))
    @constraint(m, wi_lb[e = d.Branches], Iij[e] >= -d.vmax[d.fbus[e]]*d.vmax[d.tbus[e]]*sin(d.amax[e]))

    @constraint(m, pf_lb[e = d.Branches],  pf[e] >= -d.ra[e])
    @constraint(m, pt_lb[e = d.Branches],  pt[e] >= -d.ra[e])
    @constraint(m, qf_lb[e = d.Branches],  qf[e] >= -d.ra[e])
    @constraint(m, qt_lb[e = d.Branches],  qt[e] >= -d.ra[e])
    @constraint(m, pf_ub[e = d.Branches], -pf[e] >= -d.ra[e])
    @constraint(m, pt_ub[e = d.Branches], -pt[e] >= -d.ra[e])
    @constraint(m, qf_ub[e = d.Branches], -qf[e] >= -d.ra[e])
    @constraint(m, qt_ub[e = d.Branches], -qt[e] >= -d.ra[e])

    @constraint(m, p_l_lim[g = d.Generators], pg[g] >= d.pmin[g])
    @constraint(m, p_u_lim[g = d.Generators], pg[g] <= d.pmax[g])
    
    @constraint(m, q_l_lim[g = d.Generators], qg[g] >= d.qmin[g])
    @constraint(m, q_u_lim[g = d.Generators], qg[g] <= d.qmax[g])

    @constraint(m, v_l_lim[n = d.Buses], ui[n] >= d.vmin[n]^2)
    @constraint(m, v_u_lim[n = d.Buses], ui[n] <= d.vmax[n]^2)

    return m

end

[ccoffrin]

Hi @antonhinneck! Thanks for this post. I unfortunately do not have time to work on understanding the differences between your model and the one implemented in PowerModels. However, if you try posting this injury to Julia's discourse forum others may be able to help you. See https://discourse.julialang.org/c/domain/opt/13

This implementation may also be a helpful reference for you,
https://github.com/lanl-ansi/PowerModelsAnnex.jl/blob/master/src/model/opf.jl#L137-L292

Methodologically the way I would approach debugging this would be to check the two generated JuMP models side-by-side (as you have done in this post) match up the constraints 1-by-1 and find out where the differences are. If at some point it is not clear how PowerModels arrives at some constraint or parameters I can point you to some source of where it comes from.

[antonhinneck]

Thanks @ccoffrin.
The pointer to the concise soc formulation was already a big help.
Any further questions, I will direct at the Julia discourse forum.