From 15ac519d6658e21a8ec0706f0865fcd59c7fee50 Mon Sep 17 00:00:00 2001 From: fpacaud Date: Mon, 1 Mar 2021 11:59:04 +0100 Subject: [PATCH] polar: add generic jtprod evaluation - fix bug in computation of the adjoint - refactor the computation of the objective - allow to compute the Hessian of the objective using the adjoint of active_power_constraints - remove testing of jtprod on GPU before and wait for a valid implementation on the GPU --- src/Polar/Constraints/active_power.jl | 115 ++++++----- src/Polar/Constraints/constraints.jl | 82 ++++---- src/Polar/Constraints/line_flow.jl | 40 ---- src/Polar/Constraints/power_balance.jl | 2 +- src/Polar/Constraints/reactive_power.jl | 63 ------ src/Polar/Constraints/voltage_magnitude.jl | 52 +++-- src/Polar/adjoints.jl | 212 --------------------- src/Polar/cost.jl | 102 ++++++++++ src/Polar/kernels.jl | 59 +++--- src/Polar/polar.jl | 16 +- test/Evaluators/reduced_evaluator.jl | 16 +- test/Polar/autodiff.jl | 24 ++- test/Polar/gradient.jl | 14 +- test/Polar/hessian.jl | 36 ++-- 14 files changed, 324 insertions(+), 509 deletions(-) delete mode 100644 src/Polar/adjoints.jl create mode 100644 src/Polar/cost.jl diff --git a/src/Polar/Constraints/active_power.jl b/src/Polar/Constraints/active_power.jl index d0964144..c0bcce28 100644 --- a/src/Polar/Constraints/active_power.jl +++ b/src/Polar/Constraints/active_power.jl @@ -22,44 +22,46 @@ function bounds(polar::PolarForm{T, IT, VT, MT}, ::typeof(active_power_constrain return convert(VT, pq_min), convert(VT, pq_max) end -# Jacobian -function jacobian( +function _put_active_power_injection!(fr, v_m, v_a, adj_v_m, adj_v_a, adj_P, ybus_re, ybus_im) + @inbounds for c in ybus_re.colptr[fr]:ybus_re.colptr[fr+1]-1 + to = ybus_re.rowval[c] + aij = v_a[fr] - v_a[to] + cθ = ybus_re.nzval[c]*cos(aij) + sθ = ybus_im.nzval[c]*sin(aij) + adj_v_m[fr] += v_m[to] * (cθ + sθ) * adj_P + adj_v_m[to] += v_m[fr] * (cθ + sθ) * adj_P + + adj_aij = -(v_m[fr]*v_m[to]*(ybus_re.nzval[c]*sin(aij))) + adj_aij += v_m[fr]*v_m[to]*(ybus_im.nzval[c]*cos(aij)) + adj_aij *= adj_P + adj_v_a[to] += -adj_aij + adj_v_a[fr] += adj_aij + end +end + +# Adjoint +function adjoint!( polar::PolarForm, ::typeof(active_power_constraints), - i_cons, - ∂jac, - buffer + pg, ∂pg, + vm, ∂vm, + va, ∂va, + pinj, ∂pinj, + qinj, ∂qinj, ) nbus = PS.get(polar.network, PS.NumberOfBuses()) nref = PS.get(polar.network, PS.NumberOfSlackBuses()) - index_pv = polar.indexing.index_pv - index_pq = polar.indexing.index_pq index_ref = polar.indexing.index_ref - pv_to_gen = polar.indexing.index_pv_to_gen + ref_to_gen = polar.indexing.index_ref_to_gen ybus_re, ybus_im = get(polar.topology, PS.BusAdmittanceMatrix()) - - vmag = buffer.vmag - vang = buffer.vang - adj_x = ∂jac.∇fₓ - adj_u = ∂jac.∇fᵤ - adj_vmag = ∂jac.∂vm - adj_vang = ∂jac.∂va - adj_pg = ∂jac.∂pg - fill!(adj_pg, 0.0) - fill!(adj_vmag, 0.0) - fill!(adj_vang, 0.0) - fill!(adj_x, 0.0) - fill!(adj_u, 0.0) - - adj_inj = 1.0 # Constraint on P_ref (generator) (P_inj = P_g - P_load) - bus = index_ref[i_cons] - put_active_power_injection!(bus, vmag, vang, adj_vmag, adj_vang, adj_inj, ybus_re, ybus_im) - ev = put_adjoint_kernel!(polar.device)(adj_u, adj_x, adj_vmag, adj_vang, adj_pg, - index_pv, index_pq, index_ref, pv_to_gen, - ndrange=nbus) - wait(ev) + for i in 1:nref + ibus = index_ref[i] + igen = ref_to_gen[i] + _put_active_power_injection!(ibus, vm, va, ∂vm, ∂va, ∂pg[igen], ybus_re, ybus_im) + end + return end function jacobian(polar::PolarForm, cons::typeof(active_power_constraints), buffer) @@ -86,29 +88,6 @@ function jacobian(polar::PolarForm, cons::typeof(active_power_constraints), buff return FullSpaceJacobian(jx, ju) end -# Jacobian-transpose vector product -function jtprod( - polar::PolarForm, - ::typeof(active_power_constraints), - ∂jac, - buffer, - v::AbstractVector, -) - m = size_constraint(polar, active_power_constraints) - jvx = similar(∂jac.∇fₓ) ; fill!(jvx, 0) - jvu = similar(∂jac.∇fᵤ) ; fill!(jvu, 0) - for i_cons in 1:m - if !iszero(v[i_cons]) - jacobian(polar, active_power_constraints, i_cons, ∂jac, buffer) - jx, ju = ∂jac.∇fₓ, ∂jac.∇fᵤ - jvx .+= jx .* v[i_cons] - jvu .+= ju .* v[i_cons] - end - end - ∂jac.∇fₓ .= jvx - ∂jac.∇fᵤ .= jvu -end - function hessian(polar::PolarForm, ::typeof(active_power_constraints), buffer, λ) ref = polar.indexing.index_ref pv = polar.indexing.index_pv @@ -129,3 +108,35 @@ function hessian(polar::PolarForm, ::typeof(active_power_constraints), buffer, λₚ .* ∂₂P.uu, ) end + +# MATPOWER Jacobian +function matpower_jacobian(polar::PolarForm, X::Union{State,Control}, ::typeof(active_power_constraints), V) + nbus = get(polar, PS.NumberOfBuses()) + pf = polar.network + ref = pf.ref ; nref = length(ref) + pv = pf.pv ; npv = length(pv) + pq = pf.pq ; npq = length(pq) + gen2bus = polar.indexing.index_generators + ngen = length(gen2bus) + Ybus = pf.Ybus + + dSbus_dVm, dSbus_dVa = _matpower_residual_jacobian(V, Ybus) + # w.r.t. state + if isa(X, State) + j11 = real(dSbus_dVa[ref, [pv; pq]]) + j12 = real(dSbus_dVm[ref, pq]) + return [ + j11 j12 ; + spzeros(length(pv), length(pv) + 2 * length(pq)) + ]::SparseMatrixCSC{Float64, Int} + # w.r.t. control + elseif isa(X, Control) + j11 = real(dSbus_dVm[ref, [ref; pv]]) + j12 = sparse(I, npv, npv) + return [ + j11 spzeros(length(ref), npv) ; + spzeros(npv, ngen) j12 + ]::SparseMatrixCSC{Float64, Int} + end +end + diff --git a/src/Polar/Constraints/constraints.jl b/src/Polar/Constraints/constraints.jl index d00fe5af..f05a1a21 100644 --- a/src/Polar/Constraints/constraints.jl +++ b/src/Polar/Constraints/constraints.jl @@ -15,6 +15,7 @@ include("reactive_power.jl") include("line_flow.jl") # Generic functions +## Adjoint function adjoint!( polar::PolarForm, func::Function, @@ -22,13 +23,13 @@ function adjoint!( stack, buffer, ) @assert is_constraint(func) - ∂pinj = similar(buffer.vmag) ; fill(∂pinj, 0) + ∂pinj = similar(buffer.vmag) ; fill!(∂pinj, 0.0) ∂qinj = nothing # TODO fill!(stack.∂vm, 0) fill!(stack.∂va, 0) adjoint!( polar, func, - ∂cons, cons, + cons, ∂cons, buffer.vmag, stack.∂vm, buffer.vang, stack.∂va, buffer.pinj, ∂pinj, @@ -36,6 +37,43 @@ function adjoint!( ) end +## Jacobian-transpose vector product +function jtprod( + polar::PolarForm, + func::Function, + stack, + buffer, + v::AbstractVector, +) + @assert is_constraint(func) + + m = size_constraint(polar, func) + # Adjoint w.r.t. vm, va, pinj, qinj + ∂pinj = similar(buffer.vmag) ; fill!(∂pinj, 0.0) + ∂qinj = nothing # TODO + fill!(stack.∂vm, 0) + fill!(stack.∂va, 0) + cons = similar(buffer.vmag, m) ; fill!(cons, 0.0) + adjoint!( + polar, func, + cons, v, + buffer.vmag, stack.∂vm, + buffer.vang, stack.∂va, + buffer.pinj, ∂pinj, + buffer.qinj, ∂qinj, + ) + + # Adjoint w.r.t. x and u + ∂x = stack.∇fₓ ; fill!(∂x, 0.0) + ∂u = stack.∇fᵤ ; fill!(∂u, 0.0) + adjoint_transfer!( + polar, + ∂u, ∂x, + stack.∂vm, stack.∂va, ∂pinj, ∂qinj, + ) +end + +## Sparsity detection function jacobian_sparsity(polar::PolarForm, func::Function, xx::AbstractVariable) nbus = get(polar, PS.NumberOfBuses()) Vre = Float64[i for i in 1:nbus] @@ -44,6 +82,11 @@ function jacobian_sparsity(polar::PolarForm, func::Function, xx::AbstractVariabl return matpower_jacobian(polar, xx, func, V) end +function matpower_jacobian(polar::PolarForm, func::Function, X::AbstractVariable, buffer::PolarNetworkState) + V = buffer.vmag .* exp.(im .* buffer.vang) + return matpower_jacobian(polar, X, func, V) +end + # MATPOWER's Jacobians Base.@deprecate(residual_jacobian, matpower_jacobian) @@ -78,38 +121,3 @@ function residual_jacobian(::Control, V, Ybus, pv, pq, ref) J = [j11; j21] end - -# Jacobian wrt active power generation -function active_power_jacobian(::State, V, Ybus, pv, pq, ref) - dSbus_dVm, dSbus_dVa = _matpower_residual_jacobian(V, Ybus) - j11 = real(dSbus_dVa[ref, [pv; pq]]) - j12 = real(dSbus_dVm[ref, pq]) - J = [ - j11 j12 - spzeros(length(pv), length(pv) + 2 * length(pq)) - ] -end - -function active_power_jacobian(::Control, V, Ybus, pv, pq, ref) - ngen = length(pv) + length(ref) - npv = length(pv) - dSbus_dVm, _ = _matpower_residual_jacobian(V, Ybus) - j11 = real(dSbus_dVm[ref, [ref; pv]]) - j12 = sparse(I, npv, npv) - return [ - j11 spzeros(length(ref), npv) - spzeros(npv, ngen) j12 - ] -end - -function active_power_jacobian( - polar::PolarForm, - r::AbstractVariable, - buffer::PolarNetworkState, -) - ref = polar.indexing.index_ref - pv = polar.indexing.index_pv - pq = polar.indexing.index_pq - V = buffer.vmag .* exp.(im .* buffer.vang) - return active_power_jacobian(r, V, polar.network.Ybus, pv, pq, ref) -end diff --git a/src/Polar/Constraints/line_flow.jl b/src/Polar/Constraints/line_flow.jl index 02e7b53c..807735e7 100644 --- a/src/Polar/Constraints/line_flow.jl +++ b/src/Polar/Constraints/line_flow.jl @@ -74,46 +74,6 @@ function adjoint!( wait(ev) end -# Jacobian-transpose vector product -function jtprod( - polar::PolarForm, - ::typeof(flow_constraints), - ∂jac, - buffer, - v::AbstractVector, -) - nbus = PS.get(polar.network, PS.NumberOfBuses()) - index_pv = polar.indexing.index_pv - index_pq = polar.indexing.index_pq - index_ref = polar.indexing.index_ref - pv_to_gen = polar.indexing.index_pv_to_gen - # Init buffer - adj_x = ∂jac.∇fₓ - adj_u = ∂jac.∇fᵤ - adj_vmag = ∂jac.∂vm - adj_vang = ∂jac.∂va - adj_pg = ∂jac.∂pg - ∇Jv = ∂jac.∂flow - fill!(adj_pg, 0.0) - fill!(adj_x, 0.0) - fill!(adj_u, 0.0) - fill!(∇Jv, 0.0) - # Compute gradient w.r.t. vmag and vang - flow_constraints_grad!(polar, ∇Jv, buffer, v) - # Copy results into buffer - copyto!(adj_vmag, ∇Jv[1:nbus]) - copyto!(adj_vang, ∇Jv[nbus+1:2*nbus]) - if isa(adj_x, Array) - kernel! = put_adjoint_kernel!(KA.CPU()) - else - kernel! = put_adjoint_kernel!(KA.CUDADevice()) - end - ev = kernel!(adj_u, adj_x, adj_vmag, adj_vang, adj_pg, - index_pv, index_pq, index_ref, pv_to_gen, - ndrange=nbus) - wait(ev) -end - function matpower_jacobian(polar::PolarForm, ::State, ::typeof(flow_constraints), V) nbus = get(polar, PS.NumberOfBuses()) pf = polar.network diff --git a/src/Polar/Constraints/power_balance.jl b/src/Polar/Constraints/power_balance.jl index 680bf74d..03cd81d2 100644 --- a/src/Polar/Constraints/power_balance.jl +++ b/src/Polar/Constraints/power_balance.jl @@ -104,6 +104,6 @@ function matpower_jacobian(polar::PolarForm, ::Control, ::typeof(power_balance), j12 = sparse(I, npv + npq, npv) j21 = imag(dSbus_dVm[pq, [ref; pv]]) j22 = spzeros(npq, npv) - return [j11 j12; j21 j22] + return [j11 -j12; j21 j22] end diff --git a/src/Polar/Constraints/reactive_power.jl b/src/Polar/Constraints/reactive_power.jl index 13196744..c7414898 100644 --- a/src/Polar/Constraints/reactive_power.jl +++ b/src/Polar/Constraints/reactive_power.jl @@ -87,46 +87,6 @@ function adjoint!( end # Jacobian -function jacobian( - polar::PolarForm, - ::typeof(reactive_power_constraints), - i_cons, - ∂jac, - buffer -) - nbus = PS.get(polar.network, PS.NumberOfBuses()) - gen2bus = polar.indexing.index_generators - index_pv = polar.indexing.index_pv - index_pq = polar.indexing.index_pq - index_ref = polar.indexing.index_ref - pv_to_gen = polar.indexing.index_pv_to_gen - - ybus_re, ybus_im = get(polar.topology, PS.BusAdmittanceMatrix()) - - vmag = buffer.vmag - vang = buffer.vang - adj_x = ∂jac.∇fₓ - adj_u = ∂jac.∇fᵤ - adj_vmag = ∂jac.∂vm - adj_vang = ∂jac.∂va - adj_pg = ∂jac.∂pg - fill!(adj_pg, 0.0) - fill!(adj_vmag, 0.0) - fill!(adj_vang, 0.0) - fill!(adj_x, 0.0) - fill!(adj_u, 0.0) - - adj_inj = 1.0 - bus = gen2bus[i_cons] - put_reactive_power_injection!(bus, vmag, vang, adj_vmag, adj_vang, adj_inj, ybus_re, ybus_im) - ev = put_adjoint_kernel!(polar.device)( - adj_u, adj_x, adj_vmag, adj_vang, adj_pg, - index_pv, index_pq, index_ref, pv_to_gen, - ndrange=nbus - ) - wait(ev) -end - function jacobian(polar::PolarForm, cons::typeof(reactive_power_constraints), buffer) ref = polar.indexing.index_ref pv = polar.indexing.index_pv @@ -152,29 +112,6 @@ function jacobian(polar::PolarForm, cons::typeof(reactive_power_constraints), bu return FullSpaceJacobian(jx, ju) end -# Jacobian-transpose vector product -function jtprod( - polar::PolarForm, - ::typeof(reactive_power_constraints), - ∂jac, - buffer, - v::AbstractVector, -) - m = size_constraint(polar, reactive_power_constraints) - jvx = similar(∂jac.∇fₓ) ; fill!(jvx, 0) - jvu = similar(∂jac.∇fᵤ) ; fill!(jvu, 0) - for i_cons in 1:m - if !iszero(v[i_cons]) - jacobian(polar, reactive_power_constraints, i_cons, ∂jac, buffer) - jx, ju = ∂jac.∇fₓ, ∂jac.∇fᵤ - jvx .+= jx .* v[i_cons] - jvu .+= ju .* v[i_cons] - end - end - ∂jac.∇fₓ .= jvx - ∂jac.∇fᵤ .= jvu -end - function hessian(polar::PolarForm, ::typeof(reactive_power_constraints), buffer, λ) ref = polar.indexing.index_ref pv = polar.indexing.index_pv diff --git a/src/Polar/Constraints/voltage_magnitude.jl b/src/Polar/Constraints/voltage_magnitude.jl index 6e6a358a..1a432520 100644 --- a/src/Polar/Constraints/voltage_magnitude.jl +++ b/src/Polar/Constraints/voltage_magnitude.jl @@ -1,6 +1,11 @@ is_constraint(::typeof(voltage_magnitude_constraints)) = true # We add constraint only on vmag_pq +function voltage_magnitude_constraints(polar::PolarForm, cons, vm, va, pinj, qinj) + index_pq = polar.indexing.index_pq + cons .= @view vm[index_pq] + return +end function voltage_magnitude_constraints(polar::PolarForm, cons, buffer) index_pq = polar.indexing.index_pq cons .= @view buffer.vmag[index_pq] @@ -19,17 +24,17 @@ function bounds(polar::PolarForm, ::typeof(voltage_magnitude_constraints)) return polar.x_min[fr_:to_], polar.x_max[fr_:to_] end -# State Jacobian: Jx_i = [0, ..., 1, ... 0] where -function jacobian(polar::PolarForm, ::typeof(voltage_magnitude_constraints), i_cons::Int, ∂jac, buffer) - npv = PS.get(polar.network, PS.NumberOfPVBuses()) - npq = PS.get(polar.network, PS.NumberOfPQBuses()) - fr_ = npq + npv - - # Adjoint / State - fill!(∂jac.∇fₓ, 0) - ∂jac.∇fₓ[fr_ + i_cons] = 1.0 - # Adjoint / Control - fill!(∂jac.∇fᵤ, 0) +function adjoint!( + polar::PolarForm, + ::typeof(voltage_magnitude_constraints), + cons, ∂cons, + vm, ∂vm, + va, ∂va, + pinj, ∂pinj, + qinj, ∂qinj, +) + index_pq = polar.indexing.index_pq + ∂vm[index_pq] .= ∂cons end function jacobian(polar::PolarForm, cons::typeof(voltage_magnitude_constraints), buffer) @@ -48,7 +53,7 @@ function jacobian(polar::PolarForm, cons::typeof(voltage_magnitude_constraints), return FullSpaceJacobian(jx, ju) end -function jtprod(polar::PolarForm, ::typeof(voltage_magnitude_constraints), ∂jac, buffer, v) +function jtprod(polar::PolarForm, ::typeof(voltage_magnitude_constraints), ∂jac, buffer, v::AbstractVector) npv = PS.get(polar.network, PS.NumberOfPVBuses()) npq = PS.get(polar.network, PS.NumberOfPQBuses()) fr_ = npq + npv + 1 @@ -69,3 +74,26 @@ function hessian(polar::PolarForm, ::typeof(voltage_magnitude_constraints), buff ) end +function matpower_jacobian( + polar::PolarForm, + X::Union{State,Control}, + cons::typeof(voltage_magnitude_constraints), + V, +) + m = size_constraint(polar, cons) + nᵤ = get(polar, NumberOfControl()) + nₓ = get(polar, NumberOfState()) + npv = PS.get(polar.network, PS.NumberOfPVBuses()) + npq = PS.get(polar.network, PS.NumberOfPQBuses()) + shift = npq + npv + + I = 1:m + J = (shift+1):(shift+npq) + V = ones(m) + if isa(X, State) + return sparse(I, J, V, m, nₓ) + elseif isa(X, Control) + return spzeros(m, nᵤ) + end +end + diff --git a/src/Polar/adjoints.jl b/src/Polar/adjoints.jl deleted file mode 100644 index 5e30ad0d..00000000 --- a/src/Polar/adjoints.jl +++ /dev/null @@ -1,212 +0,0 @@ -# Adjoints needed in polar formulation -# - -""" - AdjointStackObjective{VT} - -An object for storing the adjoint stack in the adjoint objective computation. - -""" -struct AdjointStackObjective{VT} - ∇fₓ::VT - ∇fᵤ::VT - ∂pg::VT - ∂vm::VT - ∂va::VT - jvₓ::VT - jvᵤ::VT - ∂flow::VT -end - -function put_active_power_injection!(fr, v_m, v_a, adj_v_m, adj_v_a, adj_P, ybus_re, ybus_im) - @inbounds for c in ybus_re.colptr[fr]:ybus_re.colptr[fr+1]-1 - to = ybus_re.rowval[c] - aij = v_a[fr] - v_a[to] - cθ = ybus_re.nzval[c]*cos(aij) - sθ = ybus_im.nzval[c]*sin(aij) - adj_v_m[fr] += v_m[to] * (cθ + sθ) * adj_P - adj_v_m[to] += v_m[fr] * (cθ + sθ) * adj_P - - adj_aij = -(v_m[fr]*v_m[to]*(ybus_re.nzval[c]*sin(aij))) - adj_aij += v_m[fr]*v_m[to]*(ybus_im.nzval[c]*cos(aij)) - adj_aij *= adj_P - adj_v_a[to] += -adj_aij - adj_v_a[fr] += adj_aij - end -end - -function put_reactive_power_injection!(fr, v_m, v_a, adj_v_m, adj_v_a, adj_P, ybus_re, ybus_im) - @inbounds for c in ybus_re.colptr[fr]:ybus_re.colptr[fr+1]-1 - to = ybus_re.rowval[c] - aij = v_a[fr] - v_a[to] - cθ = ybus_im.nzval[c]*cos(aij) - sθ = ybus_re.nzval[c]*sin(aij) - adj_v_m[fr] += v_m[to] * (sθ - cθ) * adj_P - adj_v_m[to] += v_m[fr] * (sθ - cθ) * adj_P - - adj_aij = v_m[fr]*v_m[to]*(ybus_re.nzval[c]*cos(aij)) - adj_aij += v_m[fr]*v_m[to]*(ybus_im.nzval[c]*(sin(aij))) - adj_aij *= adj_P - adj_v_a[to] += -adj_aij - adj_v_a[fr] += adj_aij - end -end - -function put!( - polar::PolarForm{T, IT, VT, MT}, - ::State, - ∂x::VT, - ∂vmag::VT, - ∂vang::VT, -) where {T, IT, VT, MT} - npv = PS.get(polar.network, PS.NumberOfPVBuses()) - npq = PS.get(polar.network, PS.NumberOfPQBuses()) - # build vector x - for (i, p) in enumerate(polar.network.pv) - ∂x[i] = ∂vang[p] - end - for (i, p) in enumerate(polar.network.pq) - ∂x[npv+i] = ∂vang[p] - ∂x[npv+npq+i] = ∂vmag[p] - end -end - -function put!( - polar::PolarForm{T, IT, VT, MT}, - ::Control, - ∂u::VT, - ∂vmag::VT, - ∂pbus::VT, -) where {T, IT, VT, MT} - npv = PS.get(polar.network, PS.NumberOfPVBuses()) - nref = PS.get(polar.network, PS.NumberOfSlackBuses()) - # build vector u - for (i, p) in enumerate(polar.network.ref) - ∂u[i] = ∂vmag[p] - end - for (i, p) in enumerate(polar.network.pv) - ∂u[nref + i] = ∂vmag[p] - ∂u[nref + npv + i] = 0.0 - end -end - -KA.@kernel function put_adjoint_kernel!( - adj_u, adj_x, adj_vmag, adj_vang, adj_pg, - index_pv, index_pq, index_ref, pv_to_gen, -) - i = @index(Global, Linear) - npv = length(index_pv) - npq = length(index_pq) - nref = length(index_ref) - - # PQ buses - if i <= npq - bus = index_pq[i] - adj_x[npv+i] = adj_vang[bus] - adj_x[npv+npq+i] = adj_vmag[bus] - # PV buses - elseif i <= npq + npv - i_ = i - npq - bus = index_pv[i_] - i_gen = pv_to_gen[i_] - adj_u[nref + i_] = adj_vmag[bus] - adj_u[nref + npv + i_] += adj_pg[i_gen] - adj_x[i_] = adj_vang[bus] - # SLACK buses - elseif i <= npq + npv + nref - i_ = i - npq - npv - bus = index_ref[i_] - adj_u[i_] = adj_vmag[bus] - end -end - -function put( - polar::PolarForm{T, IT, VT, MT}, - ::PS.Generators, - ::PS.ActivePower, - obj_autodiff::AdjointStackObjective, - buffer::PolarNetworkState -) where {T, IT, VT, MT} - ngen = PS.get(polar.network, PS.NumberOfGenerators()) - nbus = PS.get(polar.network, PS.NumberOfBuses()) - npv = PS.get(polar.network, PS.NumberOfPVBuses()) - npq = PS.get(polar.network, PS.NumberOfPQBuses()) - nref = PS.get(polar.network, PS.NumberOfSlackBuses()) - - ybus_re, ybus_im = get(polar.topology, PS.BusAdmittanceMatrix()) - - index_pv = polar.indexing.index_pv - index_ref = polar.indexing.index_ref - index_pq = polar.indexing.index_pq - pv_to_gen = polar.indexing.index_pv_to_gen - ref_to_gen = polar.indexing.index_ref_to_gen - - # Get voltages. This is only needed to get the size of adjvmag and adjvang - vmag = buffer.vmag - vang = buffer.vang - - adj_pg = obj_autodiff.∂pg - adj_x = obj_autodiff.∇fₓ - adj_u = obj_autodiff.∇fᵤ - adj_vmag = obj_autodiff.∂vm - adj_vang = obj_autodiff.∂va - fill!(adj_vmag, 0.0) - fill!(adj_vang, 0.0) - fill!(adj_x, 0.0) - fill!(adj_u, 0.0) - - for i in 1:nref - bus = index_ref[i] - i_gen = ref_to_gen[i] - # pg[i] = inj + polar.active_load[bus] - adj_inj = adj_pg[i_gen] - put_active_power_injection!(bus, vmag, vang, adj_vmag, adj_vang, adj_inj, ybus_re, ybus_im) - end - if isa(adj_x, Array) - kernel! = put_adjoint_kernel!(KA.CPU()) - else - kernel! = put_adjoint_kernel!(KA.CUDADevice()) - end - ev = kernel!(adj_u, adj_x, adj_vmag, adj_vang, adj_pg, - index_pv, index_pq, index_ref, pv_to_gen, - ndrange=nbus) - wait(ev) - - return -end - -function ∂cost(polar::PolarForm, ∂obj::AdjointStackObjective, buffer::PolarNetworkState) - pg = buffer.pg - coefs = polar.costs_coefficients - c3 = @view coefs[:, 3] - c4 = @view coefs[:, 4] - # Return adjoint of quadratic cost - ∂obj.∂pg .= c3 .+ 2.0 .* c4 .* pg - put(polar, PS.Generators(), PS.ActivePower(), ∂obj, buffer) - return -end - -function hessian_cost(polar::PolarForm, buffer::PolarNetworkState) - coefs = polar.costs_coefficients - c3 = @view coefs[:, 3] - c4 = @view coefs[:, 4] - # Change ordering from pg to [ref; pv] - pv2gen = polar.indexing.index_pv_to_gen - ref2gen = polar.indexing.index_ref_to_gen - ∂pg = (c3 .+ 2.0 .* c4 .* buffer.pg)[ref2gen] - ∂²pg = 2.0 .* c4[[ref2gen; pv2gen]] - - # Evaluate Hess-vec of objective function f(x, u) - ∂₂P = active_power_hessian(polar, buffer)::FullSpaceHessian{SparseMatrixCSC{Float64, Int}} - - ∂Pₓ = active_power_jacobian(polar, State(), buffer)::SparseMatrixCSC{Float64, Int} - ∂Pᵤ = active_power_jacobian(polar, Control(), buffer)::SparseMatrixCSC{Float64, Int} - - D = Diagonal(∂²pg) - ∇²fₓₓ = ∂pg .* ∂₂P.xx + ∂Pₓ' * D * ∂Pₓ - ∇²fᵤᵤ = ∂pg .* ∂₂P.uu + ∂Pᵤ' * D * ∂Pᵤ - ∇²fₓᵤ = ∂pg .* ∂₂P.xu + ∂Pᵤ' * D * ∂Pₓ - - return FullSpaceHessian(∇²fₓₓ, ∇²fₓᵤ, ∇²fᵤᵤ) -end - diff --git a/src/Polar/cost.jl b/src/Polar/cost.jl new file mode 100644 index 00000000..f57f33ec --- /dev/null +++ b/src/Polar/cost.jl @@ -0,0 +1,102 @@ +# Adjoints needed in polar formulation +# + +""" + AdjointStackObjective{VT} + +An object for storing the adjoint stack in the adjoint objective computation. + +""" +struct AdjointStackObjective{VT} + ∇fₓ::VT + ∇fᵤ::VT + ∂pg::VT + ∂vm::VT + ∂va::VT + jvₓ::VT + jvᵤ::VT + ∂flow::VT +end + +function cost_production(polar::PolarForm, pg) + ngen = PS.get(polar, PS.NumberOfGenerators()) + coefs = polar.costs_coefficients + c2 = @view coefs[:, 2] + c3 = @view coefs[:, 3] + c4 = @view coefs[:, 4] + # Return quadratic cost + # NB: this operation induces three allocations on the GPU, + # but is faster than writing the sum manually + cost = sum(c2 .+ c3 .* pg .+ c4 .* pg.^2) + return cost +end + +function put( + polar::PolarForm{T, IT, VT, MT}, + ::PS.Generators, + ::PS.ActivePower, + obj_autodiff::AdjointStackObjective, + buffer::PolarNetworkState +) where {T, IT, VT, MT} + + index_pv = polar.indexing.index_pv + pv_to_gen = polar.indexing.index_pv_to_gen + + adj_pg = obj_autodiff.∂pg + adj_x = obj_autodiff.∇fₓ + adj_u = obj_autodiff.∇fᵤ + adj_vmag = obj_autodiff.∂vm + adj_vang = obj_autodiff.∂va + # TODO + adj_pinj = similar(adj_vmag) ; fill!(adj_pinj, 0.0) + fill!(adj_vmag, 0.0) + fill!(adj_vang, 0.0) + + # Adjoint w.r.t Slack nodes + adjoint!(polar, active_power_constraints, adj_pg, buffer.pg, obj_autodiff, buffer) + # Adjoint w.t.t. PV nodes + adj_pinj[index_pv] .= adj_pg[pv_to_gen] + + # Adjoint w.r.t. x and u + fill!(adj_x, 0.0) + fill!(adj_u, 0.0) + adjoint_transfer!(polar, adj_u, adj_x, adj_vmag, adj_vang, adj_pinj, nothing) + + return +end + +function ∂cost(polar::PolarForm, ∂obj::AdjointStackObjective, buffer::PolarNetworkState) + pg = buffer.pg + coefs = polar.costs_coefficients + c3 = @view coefs[:, 3] + c4 = @view coefs[:, 4] + # Return adjoint of quadratic cost + ∂obj.∂pg .= c3 .+ 2.0 .* c4 .* pg + put(polar, PS.Generators(), PS.ActivePower(), ∂obj, buffer) + return +end + +function hessian_cost(polar::PolarForm, buffer::PolarNetworkState) + coefs = polar.costs_coefficients + c3 = @view coefs[:, 3] + c4 = @view coefs[:, 4] + # Change ordering from pg to [ref; pv] + pv2gen = polar.indexing.index_pv_to_gen + ref2gen = polar.indexing.index_ref_to_gen + ∂pg = (c3 .+ 2.0 .* c4 .* buffer.pg)[ref2gen] + ∂²pg = 2.0 .* c4[[ref2gen; pv2gen]] + + # Evaluate Hess-vec of objective function f(x, u) + ∂₂P = active_power_hessian(polar, buffer)::FullSpaceHessian{SparseMatrixCSC{Float64, Int}} + + ∂Pₓ = matpower_jacobian(polar, active_power_constraints, State(), buffer)::SparseMatrixCSC{Float64, Int} + ∂Pᵤ = matpower_jacobian(polar, active_power_constraints, Control(), buffer)::SparseMatrixCSC{Float64, Int} + + D = Diagonal(∂²pg) + ∇²fₓₓ = ∂pg .* ∂₂P.xx + ∂Pₓ' * D * ∂Pₓ + ∇²fᵤᵤ = ∂pg .* ∂₂P.uu + ∂Pᵤ' * D * ∂Pᵤ + ∇²fₓᵤ = ∂pg .* ∂₂P.xu + ∂Pᵤ' * D * ∂Pₓ + + return FullSpaceHessian(∇²fₓₓ, ∇²fₓᵤ, ∇²fᵤᵤ) +end + diff --git a/src/Polar/kernels.jl b/src/Polar/kernels.jl index ad2b1c20..331c4501 100644 --- a/src/Polar/kernels.jl +++ b/src/Polar/kernels.jl @@ -299,47 +299,48 @@ function transfer!(polar::PolarForm, buffer::PolarNetworkState, u) end KA.@kernel function adj_transfer_kernel!( - adj_vmag, adj_vang, adj_pinj, adj_qinj, adj_u, pv, pq, ref, adj_pload, adj_qload + adj_u, adj_x, adj_vmag, adj_vang, adj_pinj, adj_qinj, pv, pq, ref, ) i = @index(Global, Linear) npv = length(pv) npq = length(pq) nref = length(ref) - # PV bus - if i <= npv - bus = pv[i] - adj_u[nref + i] = adj_vmag[bus] - # Adjoint of P = Pg - Pd - adj_u[nref + npv + i] = adj_pinj[bus] - adj_pload[bus] = -adj_pinj[bus] - - # REF bus - else - i_ref = i - npv - bus = ref[i_ref] - adj_u[i_ref] = adj_vmag[bus] + # PQ buses + if i <= npq + bus = pq[i] + adj_x[npv+i] = adj_vang[bus] + adj_x[npv+npq+i] = adj_vmag[bus] + # PV buses + elseif i <= npq + npv + i_ = i - npq + bus = pv[i_] + adj_u[nref + i_] = adj_vmag[bus] + adj_u[nref + npv + i_] += adj_pinj[bus] + adj_x[i_] = adj_vang[bus] + # SLACK buses + elseif i <= npq + npv + nref + i_ = i - npq - npv + bus = ref[i_] + adj_u[i_] = adj_vmag[bus] end end # Transfer values in (x, u) to buffer -function adj_transfer!(polar::PolarForm, buffer::PolarNetworkState, adj_u) - if isa(u, CUDA.CuArray) - kernel! = adj_transfer_kernel!(KA.CUDADevice()) - else - kernel! = adj_transfer_kernel!(KA.CPU()) - end - nbus = length(buffer.vmag) +function adjoint_transfer!( + polar::PolarForm, + ∂u, ∂x, + ∂vm, ∂va, ∂pinj, ∂qinj, +) + nbus = get(polar, PS.NumberOfBuses()) pv = polar.indexing.index_pv pq = polar.indexing.index_pq ref = polar.indexing.index_ref - ybus_re, ybus_im = get(polar.topology, PS.BusAdmittanceMatrix()) - ev = kernel!( - buffer.adj_vmag, nothing, buffer.adj_pinj, nothing, - adj_u, - pv, pq, ref, - polar.adj_active_load, nothing, - ndrange=(length(pv)+length(ref)), + ev = adj_transfer_kernel!(polar.device)( + ∂u, ∂x, + ∂vm, ∂va, ∂pinj, ∂qinj, + pv, pq, ref; + ndrange=nbus, ) wait(ev) end @@ -575,7 +576,7 @@ KA.@kernel function adj_reactive_power_edge_kernel!( bus = ref[i_] i_gen = ref_to_gen[i_] end - inj = 0 + inj = 0.0 @inbounds for c in ybus_re_colptr[bus]:ybus_re_colptr[bus+1]-1 to = ybus_re_rowval[c] aij = vang[bus] - vang[to] diff --git a/src/Polar/polar.jl b/src/Polar/polar.jl index 5df992b7..5efd983a 100644 --- a/src/Polar/polar.jl +++ b/src/Polar/polar.jl @@ -33,8 +33,8 @@ include("kernels.jl") include("derivatives.jl") include("hessians.jl") include("getters.jl") -include("adjoints.jl") include("Constraints/constraints.jl") +include("cost.jl") function PolarForm(pf::PS.PowerNetwork, device::KA.Device) if isa(device, KA.CPU) @@ -377,20 +377,6 @@ function powerflow( return ConvergenceStatus(converged, iter, normF, sum(linsol_iters)) end -# TODO: write up a function more efficient on GPU -function cost_production(polar::PolarForm, pg) - ngen = PS.get(polar.network, PS.NumberOfGenerators()) - coefs = polar.costs_coefficients - c2 = @view coefs[:, 2] - c3 = @view coefs[:, 3] - c4 = @view coefs[:, 4] - # Return quadratic cost - # NB: this operation induces three allocations on the GPU, - # but is faster than writing the sum manually - cost = sum(c2 .+ c3 .* pg .+ c4 .* pg.^2) - return cost -end - # Utils function Base.show(io::IO, polar::PolarForm) # Network characteristics diff --git a/test/Evaluators/reduced_evaluator.jl b/test/Evaluators/reduced_evaluator.jl index 0263833a..f60e24bc 100644 --- a/test/Evaluators/reduced_evaluator.jl +++ b/test/Evaluators/reduced_evaluator.jl @@ -86,16 +86,16 @@ fill!(cons, 0) ExaPF.constraint!(nlp, cons, u) - ## Evaluation of the transpose-Jacobian product - v = similar(cons) ; fill!(v, 0) - fill!(g, 0) - ExaPF.jtprod!(nlp, g, u, v) - @test iszero(g) - fill!(v, 1) ; fill!(g, 0) - ExaPF.jtprod!(nlp, g, u, v) - ## Evaluation of the Jacobian (only on CPU) if isa(device, CPU) + ## Evaluation of the transpose-Jacobian product + v = similar(cons) ; fill!(v, 0) + fill!(g, 0) + ExaPF.jtprod!(nlp, g, u, v) + @test iszero(g) + fill!(v, 1) ; fill!(g, 0) + ExaPF.jtprod!(nlp, g, u, v) + jac = M{Float64, 2}(undef, m, n) ExaPF.jacobian!(nlp, jac, u) # Test transpose Jacobian vector product diff --git a/test/Polar/autodiff.jl b/test/Polar/autodiff.jl index 292d4f24..dc757823 100644 --- a/test/Polar/autodiff.jl +++ b/test/Polar/autodiff.jl @@ -20,9 +20,11 @@ # solve power flow conv = powerflow(polar, jx, cache, NewtonRaphson(tol=1e-12)) + V = cache.vmag .* exp.(im .* cache.vang) - # w.r.t. State + # Test Jacobian w.r.t. State for cons in [ + ExaPF.voltage_magnitude_constraints, ExaPF.power_balance, ExaPF.reactive_power_constraints, ExaPF.flow_constraints, @@ -31,6 +33,11 @@ xjacobianAD = ExaPF.AutoDiff.Jacobian(polar, cons, State()) # Evaluate Jacobian with AD J = AutoDiff.jacobian!(polar, xjacobianAD, cache) + # Matpower Jacobian + Jmat = ExaPF.matpower_jacobian(polar, State(), cons, V) + # Evaluate Jacobian transpose vector product + tgt = rand(m) + ExaPF.jtprod(polar, cons, ∂obj, cache, tgt) # Compare with FiniteDiff function jac_fd_x(x) @@ -44,10 +51,13 @@ x = [cache.vang[pv]; cache.vang[pq]; cache.vmag[pq]] Jd = FiniteDiff.finite_difference_jacobian(jac_fd_x, x) @test isapprox(Jd, xjacobianAD.J, rtol=1e-5) + @test isapprox(Jmat, xjacobianAD.J, rtol=1e-5) + @test isapprox(∂obj.∇fₓ, xjacobianAD.J' * tgt, rtol=1e-6) end - # w.r.t. Control + # Test Jacobian w.r.t. Control for cons in [ + # ExaPF.voltage_magnitude_constraints, # TODO: handle case where Jacobian is zero ExaPF.power_balance, ExaPF.reactive_power_constraints, ExaPF.flow_constraints, @@ -56,6 +66,11 @@ ujacobianAD = ExaPF.AutoDiff.Jacobian(polar, cons, Control()) # Evaluate Jacobian with AD J = AutoDiff.jacobian!(polar, ujacobianAD, cache) + # Matpower Jacobian + Jmat = ExaPF.matpower_jacobian(polar, Control(), cons, V) + # Evaluate Jacobian transpose vector product + adj = rand(m) + ExaPF.jtprod(polar, cons, ∂obj, cache, adj) # Compare with FiniteDiff function jac_fd_u(u) @@ -69,10 +84,13 @@ u = [cache.vmag[ref]; cache.vmag[pv]; cache.pinj[pv]] Jd = FiniteDiff.finite_difference_jacobian(jac_fd_u, u) @test isapprox(Jd, ujacobianAD.J, rtol=1e-5) + @test isapprox(Jmat, ujacobianAD.J, rtol=1e-6) + @test isapprox(∂obj.∇fᵤ, ujacobianAD.J' * adj, rtol=1e-6) end # Test adjoint for cons in [ + ExaPF.voltage_magnitude_constraints, ExaPF.power_balance, ExaPF.reactive_power_constraints, ExaPF.flow_constraints, @@ -80,7 +98,7 @@ m = ExaPF.size_constraint(polar, cons) λ = rand(m) c = zeros(m) - ExaPF.adjoint!(polar, cons, c, λ, ∂obj, cache) + ExaPF.adjoint!(polar, cons, λ, c, ∂obj, cache) function test_fd(vvm) cache.vmag .= vvm[1:nbus] cache.vang .= vvm[1+nbus:2*nbus] diff --git a/test/Polar/gradient.jl b/test/Polar/gradient.jl index 83bbe543..efb469a0 100644 --- a/test/Polar/gradient.jl +++ b/test/Polar/gradient.jl @@ -73,18 +73,6 @@ const KA = KernelAbstractions grad_fd = FiniteDiff.finite_difference_gradient(reduced_cost, uk) @test isapprox(grad_fd, grad_adjoint, rtol=1e-4) end - @testset "Reduced Jacobian" begin - for cons in [ - ExaPF.voltage_magnitude_constraints, - ExaPF.active_power_constraints, - ExaPF.reactive_power_constraints - ] - m = ExaPF.size_constraint(polar, cons) - for icons in 1:m - ExaPF.jacobian(polar, cons, icons, ∂obj, cache) - end - end - end @testset "Gradient of line-flow constraints" begin # Adjoint of flow_constraints() nbus = length(cache.vmag) @@ -119,7 +107,7 @@ const KA = KernelAbstractions ExaPF.flow_constraints_grad!(polar, gradg, cache, weights) @test isapprox(adgradg, fdgradg) # TODO: The gradient is slightly off with the handwritten adjoint - @test_broken isapprox(gradg, fdgradg) + @test_broken isapprox(gradg, fdgradg) end end end diff --git a/test/Polar/hessian.jl b/test/Polar/hessian.jl index 7d08384b..f670ba89 100644 --- a/test/Polar/hessian.jl +++ b/test/Polar/hessian.jl @@ -191,18 +191,6 @@ const PS = PowerSystem # Update variables x = [va[pv] ; va[pq] ; vm[pq]] u = [vm[ref]; vm[pv]; pg[pv2gen]] - # Reset voltage - V = vm .* exp.(im * va) - # Adjoint - coefs = polar.costs_coefficients - c3 = @view coefs[:, 3] - c4 = @view coefs[:, 4] - # Return adjoint of quadratic cost - adj_cost = c3 .+ 2.0 .* c4 .* pg - # Adjoint of reference generator - ∂c = adj_cost[ref2gen][1] - ∂²c = 2.0 * c4[ref2gen][1] - ∂²cpv = 2.0 * c4[pv2gen] H_ffd = FiniteDiff.finite_difference_hessian(cost_x, [x; u]) @@ -241,16 +229,6 @@ const PS = PowerSystem # h1 (state) : xl <= x <= xu # h2 (by-product) : yl <= y <= yu # Test sequential evaluation of Hessian - local ∂₂Q - for cons in [ - ExaPF.voltage_magnitude_constraints, - ExaPF.active_power_constraints, - ExaPF.reactive_power_constraints - ] - m = ExaPF.size_constraint(polar, cons) - λq = ones(m) - ∂₂Q = ExaPF.hessian(polar, cons, cache, λq) - end μ = rand(ngen) ∂₂Q = ExaPF.hessian(polar, ExaPF.reactive_power_constraints, cache, μ) @@ -297,8 +275,7 @@ const PS = PowerSystem ] @test isapprox(projp, H * tgt) - - # Line-flow + # Hessian w.r.t. Line-flow hess_lineflow = AutoDiff.Hessian(polar, ExaPF.flow_constraints) ncons = ExaPF.size_constraint(polar, ExaPF.flow_constraints) μ = rand(ncons) @@ -321,6 +298,17 @@ const PS = PowerSystem tgt = rand(nx + nu) projp = AutoDiff.adj_hessian_prod!(polar, hess_lineflow, cache, μ, tgt) @test isapprox(projp, H_fd * tgt, rtol=1e-5) + + # Hessian w.r.t. Active-power generation + hess_pg = AutoDiff.Hessian(polar, ExaPF.active_power_constraints) + μ = rand(1) + projp = AutoDiff.adj_hessian_prod!(polar, hess_pg, cache, μ, tgt) + ∂₂P = ExaPF.hessian(polar, ExaPF.active_power_constraints, cache, μ) + H = [ + ∂₂P.xx ∂₂P.xu' ; + ∂₂P.xu ∂₂P.uu + ] + @test isapprox(projp, (H * tgt)) end end