Implement Hessian of lineflow

src/Polar/Constraints/line_flow.jl

@@ -30,7 +30,7 @@ function flow_constraints_grad!(polar::PolarForm, cons_grad, buffer, weights)
     fill!(cons_grad, 0)
     adj_vmag = @view cons_grad[1:nbus]
     adj_vang = @view cons_grad[nbus+1:2*nbus]
-    ev = adj_branch_flow_kernel!(polar.device)(weights, buffer.vmag, adj_vmag, 
+    ev = adj_branch_flow_kernel!(polar.device)(weights, buffer.vmag, adj_vmag,
             buffer.vang, adj_vang,
             PT.yff_re, PT.yft_re, PT.ytf_re, PT.ytt_re,
             PT.yff_im, PT.yft_im, PT.ytf_im, PT.ytt_im,
@@ -49,6 +49,31 @@ function bounds(polar::PolarForm{T, IT, VT, MT}, ::typeof(flow_constraints)) whe
     return convert(VT, [f_min; f_min]), convert(VT, [f_max; f_max])
 end
 
+function adjoint!(
+    polar::PolarForm,
+    ::typeof(flow_constraints),
+    cons, ∂cons,
+    vm, ∂vm,
+    va, ∂va,
+    pinj, ∂pinj,
+    qinj, ∂qinj,
+)
+    nlines = PS.get(polar.network, PS.NumberOfLines())
+    nbus = PS.get(polar.network, PS.NumberOfBuses())
+    top = polar.topology
+
+    ev = adj_branch_flow_kernel!(polar.device)(
+        ∂cons,
+        vm, ∂vm,
+        va, ∂va,
+        top.yff_re, top.yft_re, top.ytf_re, top.ytt_re,
+        top.yff_im, top.yft_im, top.ytf_im, top.ytt_im,
+        top.f_buses, top.t_buses, nlines,
+        ndrange = nlines
+    )
+    wait(ev)
+end
+
 # Jacobian-transpose vector product
 function jtprod(
     polar::PolarForm,

src/Polar/kernels.jl

@@ -113,12 +113,14 @@ function residual_polar!(F, vm, va, pinj, qinj,
     wait(ev)
 end
 
-KA.@kernel function adj_residual_edge_kernel!(F, adj_F, vm, adj_vm, va, adj_va,
-                                  colptr, rowval,
-                                  ybus_re_nzval, ybus_im_nzval,
-                                  edge_vm_from, edge_vm_to,
-                                  edge_va_from, edge_va_to,
-                                  pinj, adj_pinj, qinj, pv, pq)
+KA.@kernel function adj_residual_edge_kernel!(
+    F, adj_F, vm, adj_vm, va, adj_va,
+    colptr, rowval,
+    ybus_re_nzval, ybus_im_nzval,
+    edge_vm_from, edge_vm_to,
+    edge_va_from, edge_va_to,
+    pinj, adj_pinj, qinj, pv, pq
+)
 
     npv = size(pv, 1)
     npq = size(pq, 1)
@@ -202,9 +204,11 @@ KA.@kernel function adj_residual_node_kernel!(F, adj_F, vm, adj_vm, va, adj_va,
     end
 end
 
-function adj_residual_polar!(F, adj_F, vm, adj_vm, va, adj_va,
-                         ybus_re, ybus_im,
-                         pinj, adj_pinj, qinj, pv, pq, nbus)
+function adj_residual_polar!(
+    F, adj_F, vm, adj_vm, va, adj_va,
+    ybus_re, ybus_im,
+    pinj, adj_pinj, qinj, pv, pq, nbus
+)
     npv = length(pv)
     npq = length(pq)
     nvbus = length(vm)
@@ -813,3 +817,4 @@ KA.@kernel function adj_branch_flow_kernel!(
     adj_vang[fr_bus] += adj_Δθ
     adj_vang[to_bus] -= adj_Δθ
 end
+

test/Polar/autodiff.jl

@@ -75,6 +75,7 @@
         for cons in [
             ExaPF.power_balance,
             ExaPF.reactive_power_constraints,
+            ExaPF.flow_constraints,
         ]
             m = ExaPF.size_constraint(polar, cons)
             λ = rand(m)

test/Polar/hessian.jl

@@ -296,6 +296,31 @@ const PS = PowerSystem
             ∂₂Q.xu ∂₂Q.uu
         ]
         @test isapprox(projp, H * tgt)
+
+
+        # Line-flow
+        hess_lineflow = AutoDiff.Hessian(polar, ExaPF.flow_constraints)
+        ncons = ExaPF.size_constraint(polar, ExaPF.flow_constraints)
+        μ = rand(ncons)
+        function flow_jac_x(z)
+            x_ = z[1:nx]
+            u_ = z[1+nx:end]
+            # Transfer control
+            ExaPF.transfer!(polar, cache, u_)
+            # Transfer state (manually)
+            cache.vang[pv] .= x_[1:npv]
+            cache.vang[pq] .= x_[npv+1:npv+npq]
+            cache.vmag[pq] .= x_[npv+npq+1:end]
+            V = cache.vmag .* exp.(im .* cache.vang)
+            Jx = ExaPF.matpower_jacobian(polar, State(), ExaPF.flow_constraints, V)
+            Ju = ExaPF.matpower_jacobian(polar, Control(), ExaPF.flow_constraints, V)
+            return [Jx Ju]' * μ
+        end
+
+        H_fd = FiniteDiff.finite_difference_jacobian(flow_jac_x, [x; u])
+        tgt = rand(nx + nu)
+        projp = AutoDiff.adj_hessian_prod!(polar, hess_lineflow, cache, μ, tgt)
+        @test isapprox(projp, H_fd * tgt, rtol=1e-5)
     end
 end