Skip to content

Commit

Permalink
Add support for SCOPF (#274)
Browse files Browse the repository at this point in the history 
* add formulation for SCOPF
* add support for contingencies
* add power-flow with softmin recourse
* add TrackingCost function
* add QuadraticCost function
* add proper mapping for state and control in recourse case
* add proper tests for PolarFormRecourse

N.B.: the adjoint function of power flow with recourse
      is broken on the GPU.
  • Loading branch information
@frapac
frapac authored Sep 27, 2023
1 parent dc102f7 commit 501fe67
Show file tree
Hide file tree
Showing 15 changed files with 1,316 additions and 59 deletions.
4 changes: 2 additions & 2 deletions docs/src/man/formulations.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -155,7 +155,7 @@ as `[vmag; vang; pgen]`:
# Define dimension of the problem
julia> nbus, ngen, nlines = 3, 2, 4;
julia> stack = ExaPF.NetworkStack(nbus, ngen, nlines, 1, Vector{Float64}, Vector{Float64})
julia> stack = ExaPF.NetworkStack(nbus, ngen, nlines, 0, Vector{Float64}, Vector{Float64})
8-elements NetworkStack{Vector{Float64}}
julia> stack.input
Expand Down Expand Up @@ -206,7 +206,7 @@ the problem):
```jldoctests
julia> nbus, ngen, nlines = 4, 3, 4;
julia> stack = ExaPF.NetworkStack(nbus, ngen, nlines, 1, Vector{Float64}, Vector{Float64});
julia> stack = ExaPF.NetworkStack(nbus, ngen, nlines, 0, Vector{Float64}, Vector{Float64});
julia> stack.input .= 1:length(stack.input); # index array input
Expand Down
206 changes: 206 additions & 0 deletions src/Polar/contingencies.jl
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,206 @@


abstract type AbstractContingency end

struct LineContingency <: AbstractContingency
line_id::Int
end

struct GeneratorContingency <: AbstractContingency
gen_id::Int
end

# Modify the definition of functions in presence of contingencies

function PowerFlowBalance(
polar::AbstractPolarFormulation{T, VI, VT, MT},
contingencies::Vector{LineContingency},
) where {T, VI, VT, MT}
SMT = default_sparse_matrix(polar.device)
k = nblocks(polar)

# Check we have enough blocks to represent each contingency + base case
@assert k == length(contingencies) + 1

pf = polar.network
ngen = pf.ngen
nbus = pf.nbus
gen = pf.gen2bus
npq, npv = length(pf.pq), length(pf.pv)

# Assemble matrices for loads
Cd_tot = spdiagm(nbus, nbus, ones(nbus)) # Identity matrix
Cdp = [Cd_tot[[pf.pv ; pf.pq], :]; spzeros(npq, nbus)] # active part
Cdq = [spzeros(npq+npv, nbus) ; Cd_tot[pf.pq, :]] # reactive part
Cdp = _blockdiag(Cdp, k)
Cdq = _blockdiag(Cdq, k)

# Assemble matrices for generators
Cg_tot = sparse(gen, 1:ngen, ones(ngen), nbus, ngen)
Cg = -[Cg_tot[pf.pv, :] ; spzeros(2*npq, ngen)]
Cg = _blockdiag(Cg, k)

# Assemble matrices for power flow
## Base case
M_b = PS.get_basis_matrix(polar.network)
M_ = [-M_b[[pf.pv; pf.pq; nbus .+ pf.pq], :]]
## Contingencies
for contingency in contingencies
M_c = PS.get_basis_matrix(polar.network; remove_line=contingency.line_id)
push!(M_, -M_c[[pf.pv; pf.pq; nbus .+ pf.pq], :])
end
M = blockdiag(M_...)

return PowerFlowBalance{VT, SMT}(M, Cg, Cdp, Cdq)
end

function PowerGenerationBounds(
polar::AbstractPolarFormulation{T, VI, VT, MT},
contingencies::Vector{LineContingency},
) where {T, VI, VT, MT}
SMT = default_sparse_matrix(polar.device)
k = nblocks(polar)

# Check we have enough blocks to represent each contingency + base case
@assert k == length(contingencies) + 1

pf = polar.network
nbus = pf.nbus
ns = length(pf.ref) + length(pf.pv)

# Assemble matrices for loads
Cd_tot = spdiagm(nbus, nbus, ones(nbus))
Cdp = [Cd_tot[pf.ref, :] ; spzeros(ns, nbus)]
Cdq = [spzeros(length(pf.ref), nbus) ; Cd_tot[[pf.ref ; pf.pv], :]]
Cdp = _blockdiag(Cdp, k)
Cdq = _blockdiag(Cdq, k)

# Assemble matrices for power flow
## Base case
M_b = PS.get_basis_matrix(pf)
M_ = [-M_b[[pf.ref; nbus .+ pf.ref; nbus .+ pf.pv], :]]
## Contingencies
for contingency in contingencies
M_c = PS.get_basis_matrix(pf; remove_line=contingency.line_id)
push!(M_, -M_c[[pf.ref; nbus .+ pf.ref; nbus .+ pf.pv], :])
end
M = blockdiag(M_...)

return PowerGenerationBounds{VT, SMT}(M, Cdp, Cdq)
end

function LineFlows(
polar::AbstractPolarFormulation{T,VI,VT,MT},
contingencies::Vector{LineContingency},
) where {T,VI,VT,MT}
SMT = default_sparse_matrix(polar.device)
nlines = get(polar, PS.NumberOfLines())
## Base case
Lfp, Lfq, Ltp, Ltq = PS.get_line_flow_matrices(polar.network)
Lfp_ = [Lfp]
Lfq_ = [Lfq]
Ltp_ = [Ltp]
Ltq_ = [Ltq]
## Contingencies
for contingency in contingencies
Lfp_c, Lfq_c, Ltp_c, Ltq_c = PS.get_line_flow_matrices(polar.network; remove_line=contingency.line_id)
push!(Lfp_, Lfp_c)
push!(Lfq_, Lfq_c)
push!(Ltp_, Ltp_c)
push!(Ltq_, Ltq_c)
end

return LineFlows{VT,SMT}(
nlines,
blockdiag(Lfp_...),
blockdiag(Lfq_...),
blockdiag(Ltp_...),
blockdiag(Ltq_...),
polar.device,
)
end

function PowerFlowRecourse(
polar::PolarFormRecourse{T, VI, VT, MT},
contingencies::Vector{LineContingency};
epsilon=1e-2,
alpha=nothing,
) where {T, VI, VT, MT}
@assert polar.ncustoms > 0
SMT = default_sparse_matrix(polar.device)
k = nblocks(polar)
@assert k == length(contingencies) + 1

pf = polar.network
ngen = pf.ngen
nbus = pf.nbus
gen = pf.gen2bus
npv = length(pf.pv)
npq = length(pf.pq)
@assert npv + npq + 1 == nbus

# Assemble matrices for loads
Cd_tot = spdiagm(nbus, nbus, ones(nbus)) # Identity matrix
Cdp = [Cd_tot[[pf.ref; pf.pv ; pf.pq], :]; spzeros(npq, nbus)]
Cdq = [spzeros(nbus, nbus) ; Cd_tot[pf.pq, :]]
Cdp = _blockdiag(Cdp, k)
Cdq = _blockdiag(Cdq, k)

# Assemble matrices for generators
Cg_tot = sparse(gen, 1:ngen, ones(ngen), nbus, ngen)
Cg = -[Cg_tot[[pf.ref; pf.pv], :] ; spzeros(2*npq, ngen)]
Cg = _blockdiag(Cg, k)

# Assemble matrices for power flow
## Base case
M_b = PS.get_basis_matrix(polar.network)
M_ = [-M_b[[pf.ref; pf.pv; pf.pq; nbus .+ pf.pq], :]]
## Contingencies
for contingency in contingencies
M_c = PS.get_basis_matrix(polar.network; remove_line=contingency.line_id)
push!(M_, -M_c[[pf.ref; pf.pv; pf.pq; nbus .+ pf.pq], :])
end
M = blockdiag(M_...)

# Response ratio (by default dispatch recourse evenly)
if isnothing(alpha)
alpha = ones(ngen) ./ ngen
end
@assert length(alpha) == ngen
# Bounds
_pgmin, _pgmax = PS.bounds(polar.network, PS.Generators(), PS.ActivePower())
pgmin = repeat(_pgmin, k)
pgmax = repeat(_pgmax, k)

return PowerFlowRecourse{VT, SMT}(M, Cg, Cdp, Cdq, pgmin, pgmax, alpha, epsilon)
end

function ReactivePowerBounds(
polar::PolarFormRecourse{T, VI, VT, MT},
contingencies::Vector{LineContingency},
) where {T, VI, VT, MT}
SMT = default_sparse_matrix(polar.device)
k = nblocks(polar)
@assert k == length(contingencies) + 1

pf = polar.network
nbus = pf.nbus
ns = length(pf.ref) + length(pf.pv)

# Assemble matrices for loads
Cd_tot = spdiagm(nbus, nbus, ones(nbus)) # Identity matrix
Cdq = Cd_tot[[pf.ref ; pf.pv], :]
Cdq = _blockdiag(Cdq, nblocks(polar))

# Assemble matrices for power flow
M_b = PS.get_basis_matrix(pf)
M_ = [-M_b[[nbus .+ pf.ref; nbus .+ pf.pv], :]]
for contingency in contingencies
M_c = PS.get_basis_matrix(pf; remove_line=contingency.line_id)
push!(M_, -M_c[[nbus .+ pf.ref; nbus .+ pf.pv], :])
end
M = blockdiag(M_...)

return ReactivePowerBounds{VT, SMT}(M, Cdq)
end

33 changes: 28 additions & 5 deletions src/Polar/first_order.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -41,6 +41,29 @@ function _jacobian_sparsity(polar::AbstractPolarFormulation, func::AutoDiff.Abst
return matpower_jacobian(polar, func, v)
end

function _jacobian_sparsity(polar::PolarFormRecourse, func::AutoDiff.AbstractExpression)
nbus = get(polar, PS.NumberOfBuses())
ngen = get(polar, PS.NumberOfGenerators())
v = PS.voltage(polar.network) .+ 0.01 .* rand(ComplexF64, nbus)
J1 = matpower_jacobian(polar, func, v)
m, n = size(J1)
index_gen = (2nbus+1):(2nbus+ngen)
if n == 2nbus + ngen + 1
return [J1 J1[:, index_gen]]
elseif n == 2nbus + ngen
return [J1 spzeros(m, 1) J1[:, index_gen]]
else
error("Jacobian has wrong input dimension")
end
end

function _jacobian_sparsity(polar::PolarFormRecourse, func::MultiExpressions)
return vcat([_jacobian_sparsity(polar, expr) for expr in func.exprs]...)
end
function _jacobian_sparsity(polar::PolarFormRecourse, func::ComposedExpressions)
return _jacobian_sparsity(polar, func.outer)
end

function _get_jacobian_colors(polar::AbstractPolarFormulation, func::AutoDiff.AbstractExpression, map::Vector{Int})
# Sparsity pattern
J = _jacobian_sparsity(polar, func)
Expand All @@ -51,7 +74,7 @@ function _get_jacobian_colors(polar::AbstractPolarFormulation, func::AutoDiff.Ab
end

function Jacobian(
polar::PolarForm{T, VI, VT, MT},
polar::AbstractPolarFormulation{T, VI, VT, MT},
func::AutoDiff.AbstractExpression,
map::Vector{Int},
) where {T, VI, VT, MT}
Expand All @@ -75,7 +98,7 @@ function Jacobian(
J = J_host |> SMT

# Structures
stack = NetworkStack(nbus, ngen, nlines, VT, VD)
stack = NetworkStack(nbus, ngen, nlines, polar.ncustoms, VT, VD)
init!(polar, stack)
t1sF = zeros(Float64, n_cons) |> VD

Expand All @@ -90,7 +113,7 @@ function Jacobian(

return jac
end
Jacobian(polar::PolarForm, func::AutoDiff.AbstractExpression, x::AbstractVariable) = Jacobian(polar, func, mapping(polar, x))
Jacobian(polar::AbstractPolarFormulation, func::AutoDiff.AbstractExpression, x::AbstractVariable) = Jacobian(polar, func, mapping(polar, x))


struct ArrowheadJacobian{Model, Func, Stack, VD, SMT, VI} <: AutoDiff.AbstractJacobian
Expand Down Expand Up @@ -129,7 +152,7 @@ function jacobian_arrowhead_sparsity(J, block_id, nx, nu, nblocks)
end

function ArrowheadJacobian(
polar::BlockPolarForm{T, VI, VT, MT},
polar::AbstractPolarFormulation{T, VI, VT, MT},
func::AutoDiff.AbstractExpression,
X::AbstractVariable,
) where {T, VI, VT, MT}
Expand Down Expand Up @@ -209,7 +232,7 @@ function ArrowheadJacobian(
coloring = repeat(coloring, k) |> VI

# Structures
stack = NetworkStack(nbus, ngen, nlines, k, VT, VD)
stack = NetworkStack(nbus, ngen, nlines, polar.ncustoms, k, VT, VD)
init!(polar, stack)
t1sF = zeros(Float64, n_cons) |> VD

Expand Down
Loading

0 comments on commit 501fe67

Please sign in to comment.