Problem in Power Flow using Julia


#1

I started writing “Power Flow” in two software, GAMS and Julialang, simultaneously. In addition to i did “Power Flow” in GAMS, i receive to results of “State Estimation” problem, but in Julia (using JuMP) yet i couldn’t get a feasible “Power Flow” solution.
It is possible for someone to help me in “Power Flow” in Julia?
Please do not introduce PowerModels package.


JuMP variable definition
#2

@a.shefaei this question needs more information. You are more likely to get a response if you provide a short code example and actual vs. expected result, or give some other information about your attempts and how they may have failed.


#3

two type of my code example in JuMP there are in JuMP variable definition post.


#4

That is my infeasible AC Power Flow

using JuMP
PF=Model()

# system MVA base
@NLparameter(PF, baseMVA == 100.0)
# bus data
#	bus_i	type	Pd	Qd	Gs	Bs	area	Vm	Va	baseKV	zone	Vmax	Vmin
bus = [
1        3        0        0        0        0        1        1        0        345        1        1.1        0.9
2        2        0        0        0        0        1        1        0        345        1        1.1        0.9
3        2        0        0        0        0        1        1        0        345        1        1.1        0.9
4        1        0        0        0        0        1        1        0        345        1        1.1        0.9
5        1        90       30       0        0        1        1        0        345        1        1.1        0.9
6        1        0        0        0        0        1        1        0        345        1        1.1        0.9
7        1        100      35       0        0        1        1        0        345        1        1.1        0.9
8        1        0        0        0        0        1        1        0        345        1        1.1        0.9
9        1        125      50       0        0        1        1        0        345        1        1.1        0.9
];
# generator data
#	bus	Pg	Qg	Qmax	Qmin	Vg	mBase	status	Pmax	Pmin	Pc1	Pc2	Qc1min	Qc1max	Qc2min	Qc2max	ramp_agc	ramp_10	ramp_30	ramp_q	apf
gen = [
1        0        0        300        -300        1        100        1        250        10        0        0        0        0        0        0        0        0        0        0        0
2        163      0        300        -300        1        100        1        300        10        0        0        0        0        0        0        0        0        0        0        0
3        85       0        300        -300        1        100        1        270        10        0        0        0        0        0        0        0        0        0        0        0;
];
# branch data
#	fbus	tbus	r	x	b	rateA	rateB	rateC	ratio	angle	status	angmin	angmax
branch = [
1       4       0        0.0576        0        250        250        250        0        0        1        -360        360
4       5       0.017    0.092         0.158    250        250        250        0        0        1        -360        360
5       6       0.039    0.17          0.358    150        150        150        0        0        1        -360        360
3       6       0        0.0586        0        300        300        300        0        0        1        -360        360
6       7       0.0119   0.1008        0.209    150        150        150        0        0        1        -360        360
7       8       0.0085   0.072         0.149    250        250        250        0        0        1        -360        360
8       2       0        0.0625        0        250        250        250        0        0        1        -360        360
8       9       0.032    0.161         0.306    250        250        250        0        0        1        -360        360
9       4       0.01     0.085         0.176    250        250        250        0        0        1        -360        360
];

nBus=size(bus,1); # Number of Busses

@variables PF begin
    V[i=1:nBus], (lowerbound = bus[:,13][i], start = bus[:,8][i], upperbound = bus[:,12][i])
    δ[i=1:nBus], (lowerbound = -pi, start = bus[:,9][i], upperbound = pi)
    P[i=1:nBus]
    Q[i=1:nBus]
end

# Reference Bus
@constraint(PF, constr[i=1:nBus; bus[:,2][i] == 3], V[i] == gen[:,6][i])
@constraint(PF, constr[i=1:nBus; bus[:,2][i] == 3], δ[i] == 0)
@constraint(PF, constr[i=1:nBus; bus[:,2][i] == 3], gen[:,10][i]/100 <= P[i] <= gen[:,9][i]/100)
@constraint(PF, constr[i=1:nBus; bus[:,2][i] == 3], gen[:,5][i]/100 <= Q[i] <= gen[:,4][i]/100)
# PV Busses
@constraint(PF, constr[i=1:nBus; bus[:,2][i] == 2], V[i] == gen[:,6][i]/100)
@constraint(PF, constr[i=1:nBus; bus[:,2][i] == 2], P[i] == gen[:,2][i]/100)
@constraint(PF, constr[i=1:nBus; bus[:,2][i] == 2], gen[:,5][i]/100 <= Q[i] <= gen[:,4][i]/100)
# PQ Busses
@constraint(PF, constr[i=1:nBus; bus[:,2][i] == 1], P[i] == -bus[:,3][i]/100)
@constraint(PF, constr[i=1:nBus; bus[:,2][i] == 1], Q[i] == -bus[:,4][i]/100)
@constraint(PF, constr[i=1:nBus; bus[:,2][i] == 1], bus[:,13][i] <= V[i] <= bus[:,12][i])

Y=[
17.361       0           0           17.361        0           0           0          0         0
0         16.000         0             0           0           0           0        16.000      0
0            0         17.065          0           0         17.065        0          0         0
17.361       0           0           39.448      10.689        0           0          0       11.684
0            0           0           10.689      16.166       5.733        0          0         0
0            0         17.065          0          5.733      32.246       9.852       0         0 
0            0           0             0           0          9.852      23.467     13.793      0
0          16.000        0             0           0           0         13.793     35.556     6.092
0            0           0           11.684        0           0           0         6.092    17.525
];
θ=[
-1.571     0            0         1.571        0           0           0          0           0
  0      -1.571         0          0           0           0           0         1.571        0
  0        0          -1.571       0           0           0          1.571       0           0
 1.571     0            0        -1.487      -1.388        0           0          0         -1.454
  0        0            0        -1.388      -1.370      -1.345        0          0           0
  0        0           1.571       0         -1.345      -1.495      -1.453       0           0
  0        0            0          0           0         -1.453      -1.452     -1.453        0
  0       1.571         0          0           0           0         -1.453     -1.492      -1.375
  0        0            0        -1.454        0           0           0        -1.375      -1.425
];

@NLconstraint(PF, myconstr[i=1:nBus], P[i]==V[i]*sum(V[j]*Y[i,j]*cos(δ[j]-δ[i]+θ[i,j]) for j=1:nBus))
@NLconstraint(PF, myconstr[i=1:nBus], Q[i]==V[i]*sum(V[j]*Y[i,j]*sin(δ[j]-δ[i]+θ[i,j]) for j=1:nBus))
status = solve(PF)

println("V = ", getvalue(V), " δ = ", getvalue(δ), " P = ", getvalue(P), " Q = ", getvalue(Q))

[edit: staff merged example code from other topic]


#5

I assume this was because you wanted to learn how to do this “from scratch”. That’s fine, but other readers may want to know more about the PowerModels package. It looks impressive.


#6

@a.shefaei, 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)

#7

thanks very much for your attention


#8

Hi,
I am trying to execute the codes of Test AC Power Flow and Test AC Optimal Power Flow from Problem in Power Flow using Julia. But it is showing that g is not defined i.e

g = branch[“g”] // g is not defiined b = branch[“b”] // b is not defiined c = branch[“br_b”] // it works
tr = branch[“tr”] // it works
ti = branch[“ti”] // ti is not defined
tm = tr^2 + ti^2

Could you please tell how to correct this. Currently i am using the from https://github.com/lanl-ansi/PowerModels.jl


#9

Please post a minimum reproducible example, and quote your code using backticks (`).


#10

The example posted above is for an older version of PowerModels, probably v0.2.3. The current version is v0.3.1. I’ll see if its easy to fix this example and post it here.


#11

This version of the OPF problem is valid for PowerModels v0.3.1.

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,b = PowerModels.calc_branch_y(branch)
    tr,ti = PowerModels.calc_branch_t(branch)
    c = branch["br_b"]
    tm = branch["tap"]^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)

#12

I am not an expert in power flow optimization but we solved a power flow problem using Convex.jl here. It would be interesting to know if Convex.jl can be used to satisfy your needs?
We also support optimization with complex variables and coefficients.


#13

Just to clarify this point. The referenced Convex.jl model solves the SDP relaxation of the AC Power Flow equations. In some cases the SDP model can find a globally optimal solution to the AC model, but this is not guaranteed. Also the scalability of SDP solvers is a significant challenge for power system applications.


#14

Yes, Convex.jl only solves problems which are DCP-compliant.


#15

For anyone who may be interested in this thread. From scratch JuMP models for AC-OPF like I have previously posted here are now available in this package,

A key advantage of the models in PowerModelsAnnex is that will be kept up to date with the latest version of PowerModels.