Unlike the ED model, the optimal power flow (OPF) takes into account network constraints. The latter can be achieved by using either a dc-based or ac-based power flow model. In the following, we formulate and solve both the dc-based OPF model for a simple 3 bus example using Julia's JuMP package.
The objective function of the OPF model can be formulated as follows: $$ \min \sum_{i \in I} c^g_{i} \cdot g_{i} + \sum_{b \in B} c^w \cdot w_b, $$ where $b$ is the index of buses. Note that as compared to the ED model, the OPF accounts for wind power injection at each bus separately.
s.t.
In [ ]:
using JuMP
using Interact
using Gadfly
using GraphViz
In the following experiments we employ the three bus example as in the ED and UC experiments.
In [ ]:
# Maximum power output of generators (card(I))
g_max = [1000 1500];
# Minimum power output of generators (card(I))
g_min = [0 300];
# Incremental cost of generators (card(I))
c_g = [50 100];
# Generator map (card(I) x card(b))
g_map = [0 1 0; 0 0 1];
# Incremental cost of wind generators
c_w = 50;
# Total demand (card(B))
d = [0 0 1500];
# Wind forecast (card(V))
w_f = [150 50];
# Generator map (card(I) x card(b))
w_map = [1 0 0; 0 1 0];
# Power flow limits (card(L))
f_max = [100 1000];
# Line map (card(L) x card(B))
f_map = [1 -1 0; 0 1 -1];
# Line impedance
x = [0.001 0.001]
To assess the impact of power flow limits, we will use the manipulator, varying the power flow in line L2, and observe the distribution of power flows among transmission lines. This graphical representation is enabled by package GraphViz.
In [ ]:
function writeDot(name, busIdx, busInj, renGen, f, t, lineFlow, lineLim, size=(11,14))
# This function generates a graph that richly expresses the RTS96 system state.
# name a name for the graph and resulting dot file
# busIdx bus names (could be text) in order of busInj
# busInj injection at each bus
# renGen renewable generation at each bus (0 for non-wind buses)
# f "from" node for each line
# t "to" node for each line
# lineFlow flow on each line
# lineLim limit for each line
# size size of graphical output
busInj = round(busInj,2)
lineFlow = round(lineFlow,2)
# Open the dot file, overwriting anything already there:
dotfile = IOBuffer()
# Begin writing the dot file:
write(dotfile, "digraph $(name) {\nrankdir=LR;\n")
# Set graph properties:
write(dotfile,
"graph [fontname=helvetica, tooltip=\" \", overlap=false, size=\"$(size[1]),$(size[2])\", ratio=fill, orientation=\"portrait\",layout=dot];\n")
# Set default node properties:
write(dotfile, "node [fontname=helvetica, shape=square, style=filled, fontsize=20, color=\"#bdbdbd\"];\n")
# Set default edge properties:
write(dotfile, "edge [fontname=helvetica, style=\"setlinewidth(5)\"];\n")
# Write bus data to dot file:
for i = 1:length(busIdx)
write(dotfile,
"$(i) [label=$(int(busIdx[i])), tooltip=\"Inj = $(busInj[i])\"") # bus label and tooltip
# Represent renewable nodes with blue circles:
if union(find(renGen),i) == find(renGen)
write(dotfile, ", shape=circle, color=\"#5677fc\"")
end
write(dotfile, "];\n")
end
# Write line data to file:
for i = 1:length(f)
normStrain = abs(lineFlow[i])/lineLim[i] # normalized strain on line i
# Use flow direction to determine arrow direction,
# label with flow,
# and color according to strain
if lineFlow[i] > 0
write(dotfile,
"$(f[i]) -> $(t[i]) [label=$(lineFlow[i])")
else
write(dotfile,
"$(t[i]) -> $(f[i]) [label=$(-lineFlow[i])")
end
write(dotfile,
", tooltip=\" \", labeltooltip=\"Flow = $(int(normStrain*100))%\", color=\"$(round((1 - normStrain)/3,3)) 1.000 0.700\"];\n")
end
write(dotfile, "}\n")
dottext = takebuf_string(dotfile)
#print(dottext)
return Graph(dottext)
end
In [ ]:
@manipulate for line2_limit = 1000:10:1500
f_max[2] = line2_limit
#function solve_dcopf (g_max, g_min, c_g, c_w, d, w_f)
#Define the optimal power flow (OPF) model
opf=Model()
# Define decision variables
@defVar(opf, g[i=1:2] >= 0 ) ; # power output of generators
@defVar(opf, w[v=1:2] >=0 ) ; # wind power injection
@defVar(opf, f[l=1:2]) ; # power flows
@defVar(opf, theta[b=1:3]) ; # bus angle
# Define the objective function
@setObjective(opf,Min,sum{c_g[i] * g[i],i=1:2} + sum{c_w * w[v],v=1:2});
for i in 1:2
@addConstraint(opf, g[i] <= g_max[i] ) #maximum
@addConstraint(opf, g[i] >= g_min[i] ) #minimum
end
# Define the constraint on the wind power injections
for v in 1:2
@addConstraint(opf, w[v] <= w_f[v]);
end
# Define the constraint on the power flows
for l in 1:2
@addConstraint(opf, f[l] <= f_max[l]); # direct flows
@addConstraint(opf, f[l] >= -f_max[l]); # reverse flows
end
# Define the power balance constraint
for b in 1:3
@addConstraint(opf, sum{g_map[i,b]* g[i],i=1:2} + sum{w_map[v,b] * w[v], v=1:2} + sum{f_map[l,b] * f[l], l=1:2}>= d[b]);
end
# Calculate f[l]
for l in 1:2
@addConstraint(opf, f[l] == 1/x[l] * sum{f_map [l,b] * theta[b], b=1:3}) ; # power flow in every line
end
# Slack bus
@addConstraint(opf, theta [b=1] == 0) ; # direct flows
# Solve statement
solve(opf)
writeDot("UC",[1,2,3],[0.0,getValue(g)[:]],[getValue(w)[:],0.0],[2,3],[1,2],getValue(f)[:],f_max,(3,3))
end
In [ ]: