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.
@ashefa 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.
two type of my code example in JuMP there are in JuMP variable definition post.
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]
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.
3 Likes
@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
thanks very much for your attention
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 - #6 by ccoffrin. 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 GitHub - lanl-ansi/PowerModels.jl: A Julia/JuMP Package for Power Network Optimization
Please post a minimum reproducible example, and quote your code using backticks (`).
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.
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)
1 Like
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.
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.
2 Likes
Yes, Convex.jl only solves problems which are DCP-compliant.
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.
2 Likes
Thank you ccoffrin i am following all your efforts and regarding OPF ,but i am new in this julia environment so please if you can send me energy storage constraints
Pg+Pes-Pl=B.theta
Pes.delt=E(t)-E(t-1)
-Pmax<Pes<Pmax
Emin<=Et<=Emax
The energy storage constraints are documented here, Storage Model · PowerModels
You are just getting started with Julia and JuMP you might also find this tutorial useful,
Sir,
I’m totally new to the Julia world. I tried your example code but I’m getting this error
eyError: key :bus not found
Stacktrace:
[1] getindex(::Dict{Symbol,Any}, ::Symbol) at .\dict.jl:478
[2] macro expansion at C:\Users\Biswajit Biswas.julia\packages\JuMP\tyMag\src\macros.jl:186 [inlined]
[3] top-level scope at .\In[12]:11
also, ‘network_data = PowerModels.parse_file(“$(Pkg.dir(“PowerModels”))/test/data/case30.m”)’ is not working for me, I had to give the full directory address as in ‘C:/Users/…’
Can you please help me out? Thank you.
Hi Biswajit and welcome to Julia!
The code posted in her is quite out of date at this point. Have a look at this file,
https://github.com/lanl-ansi/PowerModelsAnnex.jl/blob/master/src/model/ac-opf.jl
Which is maintained for the latest versions of Julia, JuMP, and PowerModels.