The source files for all examples can be found in /examples.
Optimal Power Flow
This example uses ConvexFlows to solve am optimal power flow problem.
The optimal power flow problem seeks a cost-minimizing plan to generate power satisfying demand in each region. We consider a network of $m$ transmission lines (edges) between $n$ regions (nodes). Each node has a demand $d_i$ and incurs cost $c_i(y_i) = (d_i - y_i)_+^2$ for power $y_i$.
Each transmission line (edge) has a loss that's a function of the power flow through it. the line and given by
\[ \ell(w) = \alpha \left(\log(1 + \exp(\beta w)) - \log 2\right) - 2w.\]
The output of a line is then $h(w) = w - \ell(w)$. The optimal power flow problem is then to minimize the total cost (maximize its negative) subject to the power flow constraints and net flow constraint:
\[\begin{aligned} & \text{maximize} && -\sum_{i=1}^n (d_i - y_i)_+^2 \\ & \text{subject to} && y = \sum_{i=1}^n A_i x_i \\ &&& h(-(x_i)_1) \le (x_i)_2 ~ \text{for} i = 1, \dots, m. \\ \end{aligned}\]
We will solve the same problem as in the standard example, but with additional optimizations.
using ConvexFlows
using Random, LinearAlgebra, SparseArrays
using Plots, LogExpFunctions, StatsBase
using Graphs: Graph, connected_components
import GraphPlot
import Cairo
const CF = ConvexFlows
const GP = GraphPlotGraphPlotGenerating the problem data
function build_graph(n; d=0.11, α=0.8, rseed=1)
Random.seed!(rseed)
xys = [sqrt(n)*rand(2) for _ in 1:n]
I = Vector{Int}()
sizehint!(I, n)
J = Vector{Int}()
sizehint!(J, n)
V = Vector{Int}()
sizehint!(V, n)
for i in 1:n, j in i+1:n
dist = norm(xys[i] - xys[j])
γ = α * min(1.0, d^2 / dist^2)
if rand() ≤ γ
push!(I, i)
push!(J, j)
push!(V, 1)
end
end
Adj = sparse(I, J, V, n, n)
Adj = Adj + Adj'
for i in 1:n
sum(Adj[i, :]) > 0 && continue
dists = [norm(xys[i] - xys[j]) for j in 1:n if j ≠ i]
jmin = argmin(dists)
Adj[i, jmin] = 1
Adj[jmin, i] = 1
end
ccs = connected_components(Graph(Adj))
nccs = length(ccs)
remaining = Set(1:nccs)
while nccs > 1
cc_i, cc_j = sample(collect(remaining), 2, replace=false)
i = rand(ccs[cc_i])
j = rand(ccs[cc_j])
Adj[i, j] = 1
Adj[j, i] = 1
ccs[cc_i] = union(ccs[cc_i], ccs[cc_j])
remaining = setdiff(remaining, Set([cc_j]))
nccs -= 1
end
return Adj, xys
end
Random.seed!(1)
n = 100
# Build a graph:
Adj, xys = build_graph(n)
network_fig_n100 = GP.gplot(
Graph(Adj),
[xy[1] for xy in xys],
[xy[2] for xy in xys],
NODESIZE=2.0/n,
nodefillc="black",
edgestrokec="dark gray",
)
Closed form solutions to arbitrage problem
# Edge gain function
h(w) = 3w - 16.0*(log1pexp(0.25 * w) - log(2))
# Closed for solution to the arbitrage problem, i.e. the w⋆ that solves
# h'(w) == ratio
wstar(ηrat, b) = ηrat ≥ 1.0 ? 0.0 : min(4.0 * log((3.0 - ηrat)/(1.0 + ηrat)), b)
lines = Edge[]
for i in 1:n, j in i+1:n
Adj[i, j] ≤ 0 && continue
# Random upper bound ub ∈ {1, 2, 3}
ub_i = rand((1., 2., 3.))
push!(lines, Edge((i, j); h=h, ub=ub_i, wstar = @inline w -> wstar(w, ub_i)))
push!(lines, Edge((j, i); h=h, ub=ub_i, wstar = @inline w -> wstar(w, ub_i)))
end
# Objective function
d = rand((0.5, 1., 2.), n)
obj = NonpositiveQuadratic(d)NonpositiveQuadratic{Float64}(100, [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0 … 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], [2.0, 0.5, 0.5, 1.0, 0.5, 2.0, 1.0, 2.0, 0.5, 2.0 … 0.5, 1.0, 1.0, 1.0, 1.0, 2.0, 1.0, 1.0, 0.5, 1.0])Solve the problem
prob = problem(obj=obj, edges=lines)
result = solve!(prob)--- Result ---
Status: OPTIMAL
f(x) : -89.41
∇f(x) : 8.95e-07
num iters : 36
solve time : 0.319s
Visualize the results
Some nodes with demand 0 are generator nodes, which we see have a net outflow.
inds = sortperm(d)
plt = bar(
prob.y[inds],
label="net flow",
color=:blue,
xticks=nothing
)
bar!(plt, d[inds], label="demand", alpha=0.5, color=:red)
Generic interface
We may also define the subproblems at each iteration directly
Arbitrage problem
We may also define the find_arb! function directly.
# ----- TransmissionLine edge ----
struct TransmissionLine{T} <: CF.Edge{T}
b::T
α::T
β::T
p_min::T
p_max::T
flow_max::T
Ai::Vector{Int}
end
function TransmissionLine(b::T, Ai::Vector{Int}) where T <: AbstractFloat
α, β = T(16), T(1/4)
p_max = one(T)
p_min = 3 - α * β * logistic(β * b)
flow_max = 3b - α * (log1pexp(β * b) - log(2))
return TransmissionLine(b, α, β, p_min, p_max, flow_max, Ai)
end
function gain(w::T, e::TransmissionLine{T}) where T
α, β, b = e.α, e.β, e.b
w = w > b ? b : w
return w - (α * (log1pexp(β * w) - log(2)) - 2w)
end
function CF.find_arb!(x::Vector{T}, e::TransmissionLine{T}, η::AbstractVector{T}) where T
η1, η2 = η[1], η[2]
x[1] = -clamp(1/e.β * log((3η2 - η1)/(η2 + η1)), zero(T), e.b)
x[2] = gain(-x[1], e)
return nothing
endObjective
struct QuadraticPowerCost{T} <: Objective
n::Int
d::Vector{T}
end
QuadraticPowerCost(d::Vector{T}) where T = QuadraticPowerCost(length(d), d)
function U(obj::QuadraticPowerCost{T}, y) where T
return -0.5*sum(x->abs2(max(x, zero(T))), obj.d - y)
end
function grad_U(obj::QuadraticPowerCost{T}, y) where T
return obj.d - y
end
function CF.Ubar(obj::QuadraticPowerCost{T}, ν) where T
return 0.5*sum(abs2, ν) - dot(obj.d, ν)
end
function CF.grad_Ubar!(g, obj::QuadraticPowerCost{T}, ν) where T
@. g = ν - obj.d
return nothing
end
CF.lower_limit(obj::QuadraticPowerCost{T}) where {T} = zeros(T, obj.n)
CF.upper_limit(obj::QuadraticPowerCost{T}) where {T} = convert(T, Inf) .+ zeros(T, obj.n)Construct problem
# Define objective function
Uy = QuadraticPowerCost(d)
# Build lines
lines_generic = TransmissionLine[]
for line in lines
push!(lines_generic, TransmissionLine(line.ub, [line.Ai[1], line.Ai[2]]))
endSolve with ConvexFlows.jl
s = Solver(
flow_objective=Uy,
edges=lines_generic,
n=n,
)
solve!(s, verbose=true)
sol_norm = norm(s.y - prob.y)
println("Norm of solution difference: $sol_norm")RUNNING THE L-BFGS-B CODE
* * *
Machine precision = 2.220D-16
N = 100 M = 5
At X0 0 variables are exactly at the bounds
At iterate 0 f= -8.00000D+01 |proj g|= 1.00000D+00
At iterate 1 f= -8.32435D+01 |proj g|= 1.84935D+00
At iterate 2 f= -8.54055D+01 |proj g|= 7.94697D-01
At iterate 3 f= -8.79083D+01 |proj g|= 6.85729D-01
At iterate 4 f= -8.85128D+01 |proj g|= 3.08336D+00
At iterate 5 f= -8.92214D+01 |proj g|= 7.38227D-01
At iterate 6 f= -8.93107D+01 |proj g|= 2.16023D-01
At iterate 7 f= -8.93609D+01 |proj g|= 1.74024D-01
At iterate 8 f= -8.93911D+01 |proj g|= 2.96564D-01
At iterate 9 f= -8.94028D+01 |proj g|= 9.86228D-02
At iterate 10 f= -8.94056D+01 |proj g|= 3.91810D-02
At iterate 11 f= -8.94069D+01 |proj g|= 1.58504D-02
At iterate 12 f= -8.94075D+01 |proj g|= 1.49024D-02
At iterate 13 f= -8.94075D+01 |proj g|= 3.43326D-02
At iterate 14 f= -8.94077D+01 |proj g|= 6.07548D-03
At iterate 15 f= -8.94077D+01 |proj g|= 3.18647D-03
At iterate 16 f= -8.94077D+01 |proj g|= 2.16905D-03
At iterate 17 f= -8.94077D+01 |proj g|= 1.04428D-02
At iterate 18 f= -8.94077D+01 |proj g|= 9.08298D-04
At iterate 19 f= -8.94077D+01 |proj g|= 5.39211D-04
At iterate 20 f= -8.94077D+01 |proj g|= 5.37826D-04
At iterate 21 f= -8.94077D+01 |proj g|= 4.63654D-04
At iterate 22 f= -8.94077D+01 |proj g|= 1.70010D-04
At iterate 23 f= -8.94077D+01 |proj g|= 7.91759D-05
At iterate 24 f= -8.94077D+01 |proj g|= 5.79295D-05
At iterate 25 f= -8.94077D+01 |proj g|= 1.78119D-04
At iterate 26 f= -8.94077D+01 |proj g|= 3.00730D-05
At iterate 27 f= -8.94077D+01 |proj g|= 1.63803D-05
At iterate 28 f= -8.94077D+01 |proj g|= 1.85695D-05
At iterate 29 f= -8.94077D+01 |proj g|= 3.67022D-05
At iterate 30 f= -8.94077D+01 |proj g|= 3.13033D-05
At iterate 31 f= -8.94077D+01 |proj g|= 3.12928D-06
* * *
Tit = total number of iterations
Tnf = total number of function evaluations
Tnint = total number of segments explored during Cauchy searches
Skip = number of BFGS updates skipped
Nact = number of active bounds at final generalized Cauchy point
Projg = norm of the final projected gradient
F = final function value
* * *
N Tit Tnf Tnint Skip Nact Projg F
100 31 33 31 0 0 3.129D-06 -8.941D+01
F = -89.407705962534408
CONVERGENCE: NORM_OF_PROJECTED_GRADIENT_<=_PGTOL
Cauchy time 1.125E-04 seconds.
Subspace minimization time 1.602E-04 seconds.
Line search time 4.226E-02 seconds.
Total User time 1.084E-01 seconds.
Norm of solution difference: 6.726324370005208e-6This page was generated using Literate.jl.