Problem in Power Flow using Julia

Specific Domains Optimization (Mathematical)
6

@ashefa, I made a quick attempt to unwrap the models that PowerModels.jl builds for AC-PF and AC-OPF. I did not test these much, but they do converge on the specified 30 bus model.

Maybe these can help you debug your code.

PS: these require the current master branch version of PowerModels.jl. But should work fine once the next version is released.

Test AC Power Flow

sing JuMP
using PowerModels
using Ipopt

ipopt_solver = IpoptSolver()

model = Model(solver = ipopt_solver)

network_data = PowerModels.parse_file("$(Pkg.dir("PowerModels"))/test/data/case30.m")
ref = PowerModels.build_ref(network_data)

@variable(model, t[i in keys(ref[:bus])])
@variable(model, v[i in keys(ref[:bus])], start = 1.0)

@variable(model, p[(l,i,j) in ref[:arcs]])
@variable(model, q[(l,i,j) in ref[:arcs]])

@variable(model, pg[i in keys(ref[:gen])])
@variable(model, qg[i in keys(ref[:gen])])


ref_bus = ref[:bus][ref[:ref_bus]]

@constraint(model, v[ref[:ref_bus]] == ref_bus["vm"])

for (i,bus) in ref[:bus]
    bus_arcs = ref[:bus_arcs][i]
    bus_gens = ref[:bus_gens][i]

    @constraint(model, sum(p[a] for a in bus_arcs) == sum{pg[g], g in bus_gens} - bus["pd"] - bus["gs"]*v[i])
    @constraint(model, sum(q[a] for a in bus_arcs) == sum{qg[g], g in bus_gens} - bus["qd"] + bus["bs"]*v[i])

    # PV Bus Constraints
    if length(ref[:bus_gens][i]) > 0 && i != ref[:ref_bus]
        # this assumes inactive generators are filtered out of bus_gens
        @assert bus["bus_type"] == 2

        @constraint(model, v[i] == bus["vm"])
        for j in ref[:bus_gens][i]
            gen = ref[:gen][j]
            @constraint(model, pg[j] == gen["pg"])
        end
    end
end

for (i,branch) in ref[:branch]
    f_bus = branch["f_bus"]
    t_bus = branch["t_bus"]
    f_idx = (i, f_bus, t_bus)
    t_idx = (i, t_bus, f_bus)

    p_fr = p[f_idx]
    p_to = p[t_idx]
    q_fr = q[f_idx]
    q_to = q[t_idx]

    v_fr = v[f_bus]
    v_to = v[t_bus]
    t_fr = t[f_bus]
    t_to = t[t_bus]

    g = branch["g"]
    b = branch["b"]
    c = branch["br_b"]
    tr = branch["tr"]
    ti = branch["ti"]
    tm = tr^2 + ti^2 

    @NLconstraint(model, p_fr == g/tm*v_fr^2 + (-g*tr+b*ti)/tm*(v_fr*v_to*cos(t_fr-t_to)) + (-b*tr-g*ti)/tm*(v_fr*v_to*sin(t_fr-t_to)) )
    @NLconstraint(model, p_to ==    g*v_to^2 + (-g*tr-b*ti)/tm*(v_to*v_fr*cos(t_to-t_fr)) + (-b*tr+g*ti)/tm*(v_to*v_fr*sin(t_to-t_fr)) )
    @NLconstraint(model, q_fr == -(b+c/2)/tm*v_fr^2 - (-b*tr-g*ti)/tm*(v_fr*v_to*cos(t_fr-t_to)) + (-g*tr+b*ti)/tm*(v_fr*v_to*sin(t_fr-t_to)) )
    @NLconstraint(model, q_to ==    -(b+c/2)*v_to^2 - (-b*tr+g*ti)/tm*(v_to*v_fr*cos(t_fr-t_to)) + (-g*tr-b*ti)/tm*(v_to*v_fr*sin(t_to-t_fr)) )
end

JuMP.solve(model)

Test AC Optimal Power Flow

using JuMP
using PowerModels
using Ipopt

ipopt_solver = IpoptSolver()

model = Model(solver = ipopt_solver)

network_data = PowerModels.parse_file("$(Pkg.dir("PowerModels"))/test/data/case30.m")
ref = PowerModels.build_ref(network_data)

@variable(model, t[i in keys(ref[:bus])])
@variable(model, v[i in keys(ref[:bus])] >= 0, start = 1.0)

@variable(model, -ref[:branch][l]["rate_a"] <= p[(l,i,j) in ref[:arcs]] <= ref[:branch][l]["rate_a"])
@variable(model, -ref[:branch][l]["rate_a"] <= q[(l,i,j) in ref[:arcs]] <= ref[:branch][l]["rate_a"])

@variable(model, ref[:gen][i]["pmin"] <= pg[i in keys(ref[:gen])] <= ref[:gen][i]["pmax"])
@variable(model, ref[:gen][i]["qmin"] <= qg[i in keys(ref[:gen])] <= ref[:gen][i]["qmax"])

@objective(model, Min, sum(gen["cost"][1]*pg[i]^2 + gen["cost"][2]*pg[i] + gen["cost"][3] for (i,gen) in ref[:gen]) )

for (i,bus) in ref[:bus]
    bus_arcs = ref[:bus_arcs][i]
    bus_gens = ref[:bus_gens][i]

    @constraint(model, sum(p[a] for a in bus_arcs) == sum{pg[g], g in bus_gens} - bus["pd"] - bus["gs"]*v[i])
    @constraint(model, sum(q[a] for a in bus_arcs) == sum{qg[g], g in bus_gens} - bus["qd"] + bus["bs"]*v[i])
end

for (i,branch) in ref[:branch]
    f_bus = branch["f_bus"]
    t_bus = branch["t_bus"]
    f_idx = (i, f_bus, t_bus)
    t_idx = (i, t_bus, f_bus)

    p_fr = p[f_idx]
    p_to = p[t_idx]
    q_fr = q[f_idx]
    q_to = q[t_idx]

    v_fr = v[f_bus]
    v_to = v[t_bus]
    t_fr = t[f_bus]
    t_to = t[t_bus]

    g = branch["g"]
    b = branch["b"]
    c = branch["br_b"]
    tr = branch["tr"]
    ti = branch["ti"]
    tm = tr^2 + ti^2 

    @NLconstraint(model, p_fr == g/tm*v_fr^2 + (-g*tr+b*ti)/tm*(v_fr*v_to*cos(t_fr-t_to)) + (-b*tr-g*ti)/tm*(v_fr*v_to*sin(t_fr-t_to)) )
    @NLconstraint(model, p_to ==    g*v_to^2 + (-g*tr-b*ti)/tm*(v_to*v_fr*cos(t_to-t_fr)) + (-b*tr+g*ti)/tm*(v_to*v_fr*sin(t_to-t_fr)) )
    @NLconstraint(model, q_fr == -(b+c/2)/tm*v_fr^2 - (-b*tr-g*ti)/tm*(v_fr*v_to*cos(t_fr-t_to)) + (-g*tr+b*ti)/tm*(v_fr*v_to*sin(t_fr-t_to)) )
    @NLconstraint(model, q_to ==    -(b+c/2)*v_to^2 - (-b*tr+g*ti)/tm*(v_to*v_fr*cos(t_fr-t_to)) + (-g*tr-b*ti)/tm*(v_to*v_fr*sin(t_to-t_fr)) )

    @constraint(model, t_fr - t_to <= branch["angmax"])
    @constraint(model, t_fr - t_to >= branch["angmin"])

    @constraint(model, p_fr^2 + q_fr^2 <= branch["rate_a"]^2)
    @constraint(model, p_to^2 + q_to^2 <= branch["rate_a"]^2)
end

JuMP.solve(model)
1 Like