diff --git a/docs/src/man/formulations.md b/docs/src/man/formulations.md
index da86bac2..c545d2e7 100644
--- a/docs/src/man/formulations.md
+++ b/docs/src/man/formulations.md
@@ -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
@@ -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
 
diff --git a/src/Polar/contingencies.jl b/src/Polar/contingencies.jl
new file mode 100644
index 00000000..39f6fcaa
--- /dev/null
+++ b/src/Polar/contingencies.jl
@@ -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
+
diff --git a/src/Polar/first_order.jl b/src/Polar/first_order.jl
index a2593573..4ed5153a 100644
--- a/src/Polar/first_order.jl
+++ b/src/Polar/first_order.jl
@@ -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)
@@ -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}
@@ -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
 
@@ -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
@@ -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}
@@ -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
 
diff --git a/src/Polar/legacy.jl b/src/Polar/legacy.jl
index 8f4bd584..5e88247d 100644
--- a/src/Polar/legacy.jl
+++ b/src/Polar/legacy.jl
@@ -29,6 +29,41 @@ function matpower_jacobian(polar::AbstractPolarFormulation, func::PowerFlowBalan
     ]::SparseMatrixCSC{Float64, Int}
 end
 
+# Warning: this method is specialized for PolarFormRecourse as it has
+# an additional state Δ
+function matpower_jacobian(polar::PolarFormRecourse, func::PowerFlowRecourse, V)
+    pf = polar.network
+    nbus = pf.nbus
+    ngen = pf.ngen
+    ref, pv, pq = pf.ref, pf.pv, pf.pq
+    gen2bus = pf.gen2bus
+    nref = length(ref)
+    npv = length(pv)
+    npq = length(pq)
+    Ybus = pf.Ybus
+
+    dSbus_dVm, dSbus_dVa = PS.matpower_residual_jacobian(V, Ybus)
+
+    Cg_tot = sparse(gen2bus, 1:ngen, ones(ngen), nbus, ngen)
+    Cg = Cg_tot[[ref; pv; pq], :]
+
+    j11 = real(dSbus_dVm[[ref; pv; pq], :])
+    j12 = real(dSbus_dVa[[ref; pv; pq], :])
+    j13 = -Cg
+    j21 = imag(dSbus_dVm[pq, :])
+    j22 = imag(dSbus_dVa[pq, :])
+    j23 = spzeros(npq, ngen)
+
+    α = func.alpha |> Array
+    j14 = -Cg * α
+    j24 = spzeros(npq, 1)
+
+    return [
+        j11 j12 j13 j14;
+        j21 j22 j23 j24
+    ]::SparseMatrixCSC{Float64, Int}
+end
+
 function matpower_jacobian(polar::AbstractPolarFormulation, func::VoltageMagnitudeBounds, V)
     pf = polar.network
     ngen = pf.ngen
@@ -89,6 +124,28 @@ function matpower_jacobian(polar::AbstractPolarFormulation, func::LineFlows, V)
     ]::SparseMatrixCSC{Float64, Int}
 end
 
+function matpower_jacobian(polar::AbstractPolarFormulation, func::ReactivePowerBounds, V)
+    pf = polar.network
+    nbus = pf.nbus
+    ngen = pf.ngen
+    gen2bus = pf.gen2bus
+    ref, pv, pq = pf.ref, pf.pv, pf.pq
+    nref = length(ref)
+    npv = length(pv)
+    npq = length(pq)
+    Ybus = pf.Ybus
+
+    dSbus_dVm, dSbus_dVa = PS.matpower_residual_jacobian(V, Ybus)
+
+    j1 = imag(dSbus_dVm[[ref; pv], :])
+    j2 = imag(dSbus_dVa[[ref; pv], :])
+    j3 = spzeros(nref + npv, ngen)
+    # w.r.t. control
+    return [
+        j1 j2 j3
+    ]::SparseMatrixCSC{Float64, Int}
+end
+
 function matpower_jacobian(polar::AbstractPolarFormulation, func::PolarBasis, V)
     pf = polar.network
     nbus = pf.nbus
@@ -185,6 +242,61 @@ function matpower_hessian(polar::AbstractPolarFormulation, func::PowerFlowBalanc
     ]::SparseMatrixCSC{Float64, Int}
 end
 
+function matpower_hessian(polar::PolarFormRecourse, func::PowerFlowRecourse, V, λ)
+    pf = polar.network
+    Ybus = pf.Ybus
+    nbus = get(polar, PS.NumberOfBuses())
+    ngen = get(polar, PS.NumberOfGenerators())
+    pq, pv, ref = pf.pq, pf.pv, pf.ref
+    npq, npv, nref = length(pq), length(pv), length(ref)
+
+    yp = zeros(nbus)
+    yp[ref] .= λ[1:nref]
+    yp[pv] .= λ[nref+1:npv+nref]
+    yp[pq] .= λ[1+npv+nref:npv+npq+nref]
+    Hpθθ, Hpvθ, Hpvv = PS._matpower_hessian(V, Ybus, yp)
+
+    yq = zeros(nbus)
+    yq[pq] .= λ[1+npv+npq:npv+2*npq]
+    Hqθθ, Hqvθ, Hqvv = PS._matpower_hessian(V, Ybus, yq)
+
+    H11 = real.(Hpvv) .+ imag.(Hqvv)
+    H12 = real.(Hpvθ) .+ imag.(Hqvθ)
+    H13 = spzeros(nbus, ngen)
+
+    H21 = real.(Hpvθ') .+ imag.(Hqvθ')
+    H22 = real.(Hpθθ) .+ imag.(Hqθθ)
+    H23 = spzeros(nbus, ngen)
+
+    H31 = spzeros(ngen, nbus)
+    H32 = spzeros(ngen, nbus)
+    H33 = spzeros(ngen, ngen)
+    return [
+        H11 H12 H13;
+        H21 H22 H23;
+        H31 H32 H33
+    ]::SparseMatrixCSC{Float64, Int}
+end
+
+function matpower_hessian(polar::PolarFormRecourse, func::QuadraticCost, V, λ)
+    pf = polar.network
+    nbus, ngen = get(polar, PS.NumberOfBuses()), get(polar, PS.NumberOfGenerators())
+    c2 = func.c2 |> Array
+    H11 = spzeros(2* nbus, 2 * nbus + ngen)
+    H21 = spzeros(ngen, 2 * nbus)
+    H22 = spdiagm(2.0 .* c2)
+    return [
+        H11;
+        H21 H22;
+    ]::SparseMatrixCSC{Float64, Int}
+end
+
+function matpower_hessian(polar::PolarFormRecourse, func::TrackingCost, V, λ)
+    pf = polar.network
+    nbus, ngen = get(polar, PS.NumberOfBuses()), get(polar, PS.NumberOfGenerators())
+    return spdiagm(ones(2nbus+ngen))
+end
+
 function matpower_hessian(polar::AbstractPolarFormulation, func::PowerGenerationBounds, V, λ)
     pf = polar.network
     Ybus = pf.Ybus
@@ -220,6 +332,37 @@ function matpower_hessian(polar::AbstractPolarFormulation, func::PowerGeneration
     ]::SparseMatrixCSC{Float64, Int}
 end
 
+function matpower_hessian(polar::AbstractPolarFormulation, func::ReactivePowerBounds, V, λ)
+    pf = polar.network
+    Ybus = pf.Ybus
+    nbus = get(polar, PS.NumberOfBuses())
+    ngen = get(polar, PS.NumberOfGenerators())
+    ref, pv = pf.ref, pf.pv
+    nref, npv = length(ref), length(pv)
+
+    yq = zeros(nbus)
+    yq[ref] .= λ[1:nref]
+    yq[pv] .= λ[nref+1:nref+npv]
+    Hqθθ, Hqvθ, Hqvv = PS._matpower_hessian(V, Ybus, yq)
+
+    H11 = imag.(Hqvv)
+    H12 = imag.(Hqvθ)
+    H13 = spzeros(nbus, ngen)
+
+    H21 = imag.(Hqvθ')
+    H22 = imag.(Hqθθ)
+    H23 = spzeros(nbus, ngen)
+
+    H31 = spzeros(ngen, nbus)
+    H32 = spzeros(ngen, nbus)
+    H33 = spzeros(ngen, ngen)
+    return [
+        H11 H12 H13;
+        H21 H22 H23;
+        H31 H32 H33
+    ]::SparseMatrixCSC{Float64, Int}
+end
+
 function matpower_hessian(polar::AbstractPolarFormulation, func::VoltageMagnitudeBounds, V, λ)
     nbus = get(polar, PS.NumberOfBuses())
     ngen = get(polar, PS.NumberOfGenerators())
diff --git a/src/Polar/newton.jl b/src/Polar/newton.jl
index edf1bce3..cd07cb2c 100644
--- a/src/Polar/newton.jl
+++ b/src/Polar/newton.jl
@@ -162,6 +162,9 @@ function nlsolve!(
 
         iter += 1
     end
+    # Final evaluation to refresh all values in stack
+    jac.func(residual, stack)
+
     time_total = time() - tic
     return ConvergenceStatus(
         converged, iter, normF, sum(linsol_iters), time_jacobian, time_linear_solver, time_total,
diff --git a/src/Polar/polar.jl b/src/Polar/polar.jl
index fe3678b5..ba1b4f7f 100644
--- a/src/Polar/polar.jl
+++ b/src/Polar/polar.jl
@@ -36,14 +36,15 @@ giving a mathematical formulation with:
 struct PolarForm{T, IT, VT, MT} <: AbstractPolarFormulation{T, IT, VT, MT}
     network::PS.PowerNetwork
     device::KA.Backend
+    ncustoms::Int  # custom variables defined by user
 end
 
-function PolarForm(pf::PS.PowerNetwork, device::KA.CPU)
-    return PolarForm{Float64, Vector{Int}, Vector{Float64}, Matrix{Float64}}(pf, device)
+function PolarForm(pf::PS.PowerNetwork, device::KA.CPU, ncustoms::Int=0)
+    return PolarForm{Float64, Vector{Int}, Vector{Float64}, Matrix{Float64}}(pf, device, ncustoms)
 end
 # Convenient constructor
-PolarForm(datafile::String, device=CPU()) = PolarForm(PS.PowerNetwork(datafile), device)
-PolarForm(polar::PolarForm, device=CPU()) = PolarForm(polar.network, device)
+PolarForm(datafile::String, device=CPU(); ncustoms=0) = PolarForm(PS.PowerNetwork(datafile), device, ncustoms)
+PolarForm(polar::PolarForm, device=CPU()) = PolarForm(polar.network, device, polar.ncustoms)
 
 introduce(polar::PolarForm) = "Polar formulation"
 nblocks(polar::PolarForm) = 1
@@ -60,12 +61,13 @@ struct BlockPolarForm{T, IT, VT, MT} <: AbstractPolarFormulation{T, IT, VT, MT}
     network::PS.PowerNetwork
     device::KA.Backend
     k::Int
+    ncustoms::Int  # custom variables defined by user
 end
-function BlockPolarForm(pf::PS.PowerNetwork, device, k::Int)
-    return BlockPolarForm{Float64, Vector{Int}, Vector{Float64}, Matrix{Float64}}(pf, device, k)
+function BlockPolarForm(pf::PS.PowerNetwork, device, k::Int, ncustoms::Int=0)
+    return BlockPolarForm{Float64, Vector{Int}, Vector{Float64}, Matrix{Float64}}(pf, device, k, ncustoms)
 end
-BlockPolarForm(datafile::String, k::Int, device=CPU()) = BlockPolarForm(PS.PowerNetwork(datafile), device, k)
-BlockPolarForm(polar::PolarForm, k::Int) = BlockPolarForm(polar.network, polar.device, k)
+BlockPolarForm(datafile::String, k::Int, device=CPU(); ncustoms=0) = BlockPolarForm(PS.PowerNetwork(datafile), device, k, ncustoms)
+BlockPolarForm(polar::PolarForm, k::Int) = BlockPolarForm(polar.network, polar.device, k, polar.ncustoms)
 
 introduce(polar::BlockPolarForm) = "$(polar.k)-BlockPolar formulation"
 nblocks(polar::BlockPolarForm) = polar.k
@@ -117,8 +119,8 @@ giving a mathematical formulation with:
 ```
 
 """
-function load_polar(case, device=CPU(); dir=PS.EXADATA)
-    return PolarForm(PS.load_case(case, dir), device)
+function load_polar(case, device=CPU(); ncustoms=0, dir=PS.EXADATA)
+    return PolarForm(PS.load_case(case, dir), device, ncustoms)
 end
 
 """
@@ -273,6 +275,8 @@ end
 
 include("stacks.jl")
 include("functions.jl")
+include("recourse.jl")
+include("contingencies.jl")
 include("first_order.jl")
 include("second_order.jl")
 include("newton.jl")
diff --git a/src/Polar/recourse.jl b/src/Polar/recourse.jl
new file mode 100644
index 00000000..e5a75b87
--- /dev/null
+++ b/src/Polar/recourse.jl
@@ -0,0 +1,543 @@
+
+struct PolarFormRecourse{T, IT, VT, MT} <: AbstractPolarFormulation{T, IT, VT, MT}
+    network::PS.PowerNetwork
+    device::KA.Backend
+    k::Int
+    ncustoms::Int  # custom variables defined by user
+end
+
+function PolarFormRecourse(pf::PS.PowerNetwork, device::KA.CPU, k::Int)
+    ngen = PS.get(pf, PS.NumberOfGenerators())
+    ncustoms = (ngen + 1) * k
+    return PolarFormRecourse{Float64, Vector{Int}, Vector{Float64}, Matrix{Float64}}(pf, device, k, ncustoms)
+end
+# Convenient constructor
+PolarFormRecourse(datafile::String, k::Int, device=CPU()) = PolarFormRecourse(PS.PowerNetwork(datafile), device, k)
+PolarFormRecourse(polar::PolarForm, k::Int) = PolarFormRecourse(polar.network, polar.device, k)
+
+name(polar::PolarFormRecourse) = "Polar formulation with recourse"
+nblocks(polar::PolarFormRecourse) = polar.k
+
+"""
+    ExtendedPowerFlowBalance
+"""
+struct PowerFlowRecourse{VT, MT} <: AutoDiff.AbstractExpression
+    M::MT
+    Cg::MT
+    Cdp::MT
+    Cdq::MT
+    pgmin::VT
+    pgmax::VT
+    alpha::VT
+    epsilon::Float64
+end
+
+function PowerFlowRecourse(
+    polar::PolarFormRecourse{T, VI, VT, MT};
+    epsilon=1e-2,
+    alpha=nothing,
+) where {T, VI, VT, MT}
+    @assert polar.ncustoms > 0
+    SMT = default_sparse_matrix(polar.device)
+    k = nblocks(polar)
+
+    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
+    Cg_tot = sparse(gen, 1:ngen, ones(ngen), nbus, ngen)
+    Cd_tot = spdiagm(nbus, nbus, ones(nbus)) # Identity matrix
+    Cg = -[Cg_tot[[pf.ref; pf.pv], :] ; spzeros(2*npq, ngen)]
+    M_tot = PS.get_basis_matrix(polar.network)
+    M = -M_tot[[pf.ref; pf.pv; pf.pq; nbus .+ pf.pq], :]
+    # constant term
+    Cdp = [Cd_tot[[pf.ref; pf.pv ; pf.pq], :]; spzeros(npq, nbus)]
+    Cdq = [spzeros(nbus, nbus) ; Cd_tot[pf.pq, :]]
+
+    M   = _blockdiag(M, k)
+    Cg  = _blockdiag(Cg, k)
+    Cdp = _blockdiag(Cdp, k)
+    Cdq = _blockdiag(Cdq, k)
+
+    # 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
+
+Base.length(func::PowerFlowRecourse) = size(func.M, 1)
+
+@inline function _softmin(x1, x2, ϵ)
+    return (x1 - ϵ * log1p(exp((x1 - x2) / ϵ)))
+end
+
+# Smooth approximation of max(pmin, min(p, pmax))
+# (numerically stable version)
+@inline function smooth_response(p, pmin, pmax, ϵ)
+    threshold = 100.0
+    if p >= pmax + threshold * ϵ
+        return pmax
+    elseif p >= 0.5 * (pmax + pmin)
+        return _softmin(p, pmax, ϵ)
+    elseif p >= (pmin - threshold * ϵ)
+        return -_softmin(-p, -pmin, ϵ)
+    else
+        return pmin
+    end
+end
+
+@kernel function smooth_pgen!(
+    pgen, @Const(setpoint), @Const(delta), @Const(alpha), @Const(pgmin), @Const(pgmax), epsilon, ngen, nblocks,
+)
+    i, j = @index(Global, NTuple)
+    idx = i + (j-1) * ngen
+    @inbounds begin
+        p = setpoint[idx] + delta[j] * alpha[i]
+        pgen[idx] = smooth_response(p, pgmin[idx], pgmax[idx], epsilon)
+    end
+end
+
+function (func::PowerFlowRecourse)(cons::AbstractArray, stack::AbstractNetworkStack)
+    k = nblocks(stack)
+    ngen = div(length(stack.pgen), k)
+    fill!(cons, 0.0)
+    # Constant terms
+    mul!(cons, func.Cdp, stack.pload, 1.0, 1.0)
+    mul!(cons, func.Cdq, stack.qload, 1.0, 1.0)
+    # Variable terms
+    mul!(cons, func.M, stack.ψ, 1.0, 1.0)
+
+    Δ = view(stack.vuser, 1:k)
+    setpoint = view(stack.vuser, k+1:k+k*ngen)
+
+    # Kernel evaluation
+    device = KA.get_backend(stack.pgen)
+    ndrange = (ngen, k)
+    smooth_pgen!(device)(
+        stack.pgen, setpoint, Δ, func.alpha, func.pgmin, func.pgmax, func.epsilon, ngen, k,
+        ndrange=ndrange,
+    )
+    KA.synchronize(device)
+
+    mul!(cons, func.Cg, stack.pgen, 1.0, 1.0)
+    return
+end
+
+# adjoint
+function adjoint_smooth_response(p::T, pmin, pmax, ϵ) where T
+    threshold = 100.0
+    if p >= pmax + threshold * ϵ
+        return 0.0
+    elseif p >= 0.5 * (pmax + pmin)
+        X = 1.0 / ϵ * (p - pmax)
+        return 1.0 / (exp(X) + 1.0)
+    elseif p >= (pmin - threshold * ϵ)
+        X = 1.0 / ϵ * (p - pmin)
+        return 1.0 / (exp(-X) + 1.0)
+    else
+        return 0.0
+    end
+end
+
+function adjoint_smooth_pgen!(
+    ∂setpoint, ∂delta,
+    setpoint, delta, pgen, ∂pgen,
+    alpha, pgmin, pgmax, epsilon, ngen, nblocks,
+)
+    for j in 1:nblocks
+        for i in 1:ngen
+            idx = i + (j-1) * ngen
+            p = setpoint[idx] + delta[j] * alpha[i]
+            ∂p = adjoint_smooth_response(p, pgmin[idx], pgmax[idx], epsilon) * ∂pgen[idx]
+            ∂delta[j] += alpha[i] * ∂p
+            ∂setpoint[idx] += ∂p
+        end
+    end
+end
+
+function adjoint!(func::PowerFlowRecourse, ∂stack, stack, ∂v)
+    k = nblocks(stack)
+    ngen = div(length(stack.pgen), k)
+    Δ = view(stack.vuser, 1:k)
+    setpoint = view(stack.vuser, k+1:k+k*ngen)
+    ∂Δ = view(∂stack.vuser, 1:k)
+    ∂setpoint = view(∂stack.vuser, k+1:k+k*ngen)
+
+    mul!(∂stack.pgen, func.Cg', ∂v, 1.0, 1.0)
+    adjoint_smooth_pgen!(
+        ∂setpoint, ∂Δ, setpoint, Δ, stack.pgen, ∂stack.pgen,
+        func.alpha, func.pgmin, func.pgmax, func.epsilon, ngen, k,
+    )
+    mul!(∂stack.ψ, func.M', ∂v, 1.0, 1.0)
+    return
+end
+
+function bounds(polar::AbstractPolarFormulation{T,VI,VT,MT}, func::PowerFlowRecourse) where {T,VI,VT,MT}
+    m = length(func)
+    return (fill!(VT(undef, m), zero(T)) , fill!(VT(undef, m), zero(T)))
+end
+
+function Base.show(io::IO, func::PowerFlowRecourse)
+    print(io, "PowerFlowRecourse (AbstractExpression)")
+end
+
+
+"""
+    ReactivePowerBounds
+"""
+struct ReactivePowerBounds{VT, MT} <: AutoDiff.AbstractExpression
+    M::MT
+    Cdq::MT
+end
+
+function ReactivePowerBounds(polar::PolarFormRecourse{T, VI, VT, MT}) where {T, VI, VT, MT}
+    SMT = default_sparse_matrix(polar.device)
+    pf = polar.network
+    nbus = pf.nbus
+    M_tot = PS.get_basis_matrix(pf)
+    ns = length(pf.ref) + length(pf.pv)
+
+    M = -M_tot[[nbus .+ pf.ref; nbus .+ pf.pv], :]
+    Cd_tot = spdiagm(nbus, nbus, ones(nbus)) # Identity matrix
+    Cdq = Cd_tot[[pf.ref ; pf.pv], :]
+
+    M = _blockdiag(M, nblocks(polar))
+    Cdq = _blockdiag(Cdq, nblocks(polar))
+
+    return ReactivePowerBounds{VT, SMT}(M, Cdq)
+end
+
+Base.length(func::ReactivePowerBounds) = size(func.M, 1)
+
+function (func::ReactivePowerBounds)(cons::AbstractArray, stack::AbstractNetworkStack)
+    fill!(cons, 0.0)
+    # Constant terms
+    mul!(cons, func.Cdq, stack.qload, 1.0, 1.0)
+    # Variable terms
+    mul!(cons, func.M, stack.ψ, 1.0, 1.0)
+    return
+end
+
+function adjoint!(func::ReactivePowerBounds, ∂stack, stack, ∂v)
+    mul!(∂stack.ψ, func.M', ∂v, 1.0, 1.0)
+    return
+end
+
+function bounds(polar::PolarFormRecourse{T,VI,VT,MT}, func::ReactivePowerBounds) where {T,VI,VT,MT}
+    pf = polar.network
+    ngen = pf.ngen
+    nbus = pf.nbus
+    ref, pv = pf.ref, pf.pv
+    # Build incidence matrix
+    Cg = sparse(pf.gen2bus, 1:ngen, ones(ngen), nbus, ngen)
+    Cgq = Cg[[ref ; pv], :]
+    # Get original bounds
+    q_min, q_max = PS.bounds(polar.network, PS.Generators(), PS.ReactivePower())
+    # Aggregate bounds on ref and pv nodes
+    lb = Cgq * q_min
+    ub = Cgq * q_max
+    return (
+        convert(VT, repeat(lb, nblocks(polar))),
+        convert(VT, repeat(ub, nblocks(polar))),
+    )
+end
+
+function Base.show(io::IO, func::ReactivePowerBounds)
+    print(io, "ReactivePowerBounds (AbstractExpression)")
+end
+
+
+#=
+    Mapping
+=#
+
+function mapping(polar::PolarFormRecourse, ::Control, k::Int=1)
+    pf = polar.network
+    nbus = get(polar, PS.NumberOfBuses())
+    ngen = get(polar, PS.NumberOfGenerators())
+    ngen = polar.network.ngen
+    ref, pv = pf.ref, pf.pv
+    genidx = 1:ngen
+    nu = (length(pv) + length(ref) + length(genidx)) * k
+    mapu = zeros(Int, nu)
+
+    shift_mag = 0
+    shift_pgen = (2 * nbus + ngen + 1) * k
+    index = 1
+    for i in 1:k
+        for j in ref
+            mapu[index] = j + (i-1) * nbus + shift_mag
+            index += 1
+        end
+        for j in pv
+            mapu[index] = j + (i-1) * nbus + shift_mag
+            index += 1
+        end
+        for j in genidx
+            mapu[index] = j + (i-1) * ngen + shift_pgen
+            index += 1
+        end
+    end
+    return mapu
+end
+
+function mapping(polar::PolarFormRecourse, ::State, k::Int=1)
+    pf = polar.network
+    nbus = get(polar, PS.NumberOfBuses())
+    ngen = get(polar, PS.NumberOfGenerators())
+    ndelta = 1
+    ref, pv, pq = pf.ref, pf.pv, pf.pq
+    nx = (length(pv) + 2*length(pq) + ndelta) * k
+    mapx = zeros(Int, nx)
+
+    shift_ang = k * nbus
+    shift_mag = 0
+    shift_delta = (2 * nbus + ngen) * k
+    index = 1
+    for i in 1:k
+        for j in [pv; pq]
+            mapx[index] = j + (i-1) * nbus + shift_ang
+            index += 1
+        end
+        for j in pq
+            mapx[index] = j + (i-1) * nbus + shift_mag
+            index += 1
+        end
+        # delta
+        mapx[index] = 1 + (i-1) * ndelta + shift_delta
+        index += 1
+    end
+    return mapx
+end
+
+function mapping(polar::PolarFormRecourse, ::AllVariables, k::Int=1)
+    pf = polar.network
+    nbus = get(polar, PS.NumberOfBuses())
+    ngen = get(polar, PS.NumberOfGenerators())
+    ndelta = 1
+    ref, pv, pq = pf.ref, pf.pv, pf.pq
+    genidx = collect(1:ngen)
+    nx = (length(pv) + 2*length(pq) + ndelta) * k
+    nu = (length(pv) + length(ref) + length(genidx)) * k
+    mapxu = zeros(Int, nx+nu)
+
+    shift_mag = 0
+    shift_ang = k * nbus
+    shift_pgen = (2 * nbus + ngen + ndelta) * k
+    shift_delta = (2 * nbus + ngen) * k
+    index = 1
+    for i in 1:k
+        # / x
+        for j in [pv; pq]
+            mapxu[index] = j + (i-1) * nbus + shift_ang
+            index += 1
+        end
+        for j in pq
+            mapxu[index] = j + (i-1) * nbus + shift_mag
+            index += 1
+        end
+        # delta
+        mapxu[index] = 1 + (i-1) * ndelta + shift_delta
+        index += 1
+        # / u
+        for j in [ref; pv]
+            mapxu[index] = j + (i-1) * nbus + shift_mag
+            index += 1
+        end
+        for j in genidx
+            mapxu[index] = j + (i-1) * ngen + shift_pgen
+            index += 1
+        end
+    end
+    return mapxu
+end
+
+#=
+    Stack
+=#
+
+function init!(polar::PolarFormRecourse, stack::NetworkStack; update_loads=true)
+    vmag = get(polar.network, PS.VoltageMagnitude())
+    vang = get(polar.network, PS.VoltageAngle())
+    pgen = get(polar.network, PS.ActivePower())
+    pload = get(polar.network, PS.ActiveLoad())
+    qload = get(polar.network, PS.ReactiveLoad())
+
+    nbus = length(vmag)
+    ngen = length(pgen)
+    nload = length(pload)
+
+    for s in 1:nblocks(stack)
+        i = (s - 1) * nbus  + 1
+        copyto!(stack.vmag, i, vmag, 1, nbus)
+        copyto!(stack.vang, i, vang, 1, nbus)
+
+        i = (s - 1) * ngen  + 1
+        copyto!(stack.pgen, i, pgen, 1, ngen)
+
+        i = (s - 1) * ngen  + nblocks(stack) + 1
+        copyto!(stack.vuser, i, pgen, 1, ngen)
+
+        if update_loads
+            i = (s - 1) * nload  + 1
+            copyto!(stack.pload, i, pload, 1, nload)
+            copyto!(stack.qload, i, qload, 1, nload)
+        end
+    end
+end
+
+function bounds(polar::PolarFormRecourse{T, VI, VT, MT}, stack::NetworkStack) where {T, VI, VT, MT}
+    nbus = polar.network.nbus
+    ngen = polar.network.ngen
+    vmag_min, vmag_max = PS.bounds(polar.network, PS.Buses(), PS.VoltageMagnitude())
+    vang_min, vang_max = fill(-Inf, nbus), fill(Inf, nbus)
+    pgen_min, pgen_max = PS.bounds(polar.network, PS.Generators(), PS.ActivePower())
+    delta_min, delta_max = -Inf, Inf
+
+    # NB: we set delta_min, delta_max = 0, 0 for the first scenario
+    # to ensure it is set to the reference value (no adjustment)
+    lb = [
+        repeat(vmag_min, nblocks(polar));
+        repeat(vang_min, nblocks(polar));
+        repeat(pgen_min, nblocks(polar));
+        [0.0; fill(delta_min, nblocks(polar)-1)];
+        pgen_min;
+    ]
+    ub = [
+        repeat(vmag_max, nblocks(polar));
+        repeat(vang_max, nblocks(polar));
+        repeat(pgen_max, nblocks(polar));
+        [0.0; fill(delta_max, nblocks(polar)-1)];
+        pgen_max;
+    ]
+    return convert(VT, lb), convert(VT, ub)
+end
+
+struct QuadraticCost{VT, MT} <: AutoDiff.AbstractExpression
+    c0::VT
+    c1::VT
+    c2::VT
+    k::Int
+    device::KA.Backend
+end
+
+function QuadraticCost(polar::PolarFormRecourse{T, VI, VT, MT}) where {T, VI, VT, MT}
+    coefs = PS.get_costs_coefficients(polar.network)
+    c0 = @view coefs[:, 2]
+    c1 = @view coefs[:, 3]
+    c2 = @view coefs[:, 4]
+    return QuadraticCost{VT, Nothing}(c0, c1, c2, nblocks(polar), polar.device)
+end
+
+Base.length(func::QuadraticCost) = func.k
+
+function (func::QuadraticCost)(output::AbstractArray, stack::AbstractNetworkStack)
+    k = nblocks(stack)
+    ngen = length(func.c0)
+    costs = stack.intermediate.c
+    setpoint = view(stack.vuser, k+1:k+k*ngen)
+    # Compute quadratic costs
+    ndrange = (ngen, k)
+    _quadratic_cost_kernel(func.device)(
+        costs, setpoint, func.c0, func.c1, func.c2, ngen;
+        ndrange=ndrange,
+    )
+    KA.synchronize(func.device)
+    output .= sum(reshape(costs, ngen, nblocks(stack))', dims=2)
+    return
+end
+
+function adjoint!(func::QuadraticCost, ∂stack, stack, ∂v)
+    k = nblocks(stack)
+    ngen = length(func.c0)
+    setpoint = view(stack.vuser, k+1:k+k*ngen)
+    ∂setpoint = view(∂stack.vuser, k+1:k+k*ngen)
+
+    ndrange = (ngen, nblocks(stack))
+    _adj_quadratic_cost_kernel(func.device)(
+        ∂setpoint, setpoint, ∂v, func.c0, func.c1, func.c2, ngen;
+        ndrange=ndrange,
+    )
+    KA.synchronize(func.device)
+    return
+end
+
+
+struct TrackingCost{VT, MT} <: AutoDiff.AbstractExpression
+    πv::Float64
+    πa::Float64
+    πp::Float64
+    vmag_ref::VT
+    vang_ref::VT
+    pgen_ref::VT
+    nbus::Int
+    ngen::Int
+    k::Int
+end
+
+function TrackingCost(
+    polar::PolarFormRecourse{T, VI, VT, MT},
+    stack_ref::AbstractNetworkStack;
+    weight_vmag=1.0, weight_vang=1.0, weight_pgen=1.0,
+)  where {T, VI, VT, MT}
+    nbus = PS.get(polar, PS.NumberOfBuses())
+    ngen = PS.get(polar, PS.NumberOfGenerators())
+    return TrackingCost{VT, MT}(
+        weight_vmag, weight_vang, weight_pgen,
+        stack_ref.vmag[1:nbus], stack_ref.vang[1:nbus], stack_ref.pgen[1:ngen],
+        nbus, ngen, nblocks(polar),
+    )
+end
+
+Base.length(func::TrackingCost) = func.k
+
+function (func::TrackingCost)(output::AbstractArray, stack::AbstractNetworkStack)
+    k = func.k
+    vm = view(stack.vmag, 1:func.nbus)
+    dvm = vm .- func.vmag_ref
+    q1 = 0.5 * func.πv * dot(dvm, dvm)
+
+    va = view(stack.vang, 1:func.nbus)
+    dva = va .- func.vang_ref
+    q2 = 0.5 * func.πa * dot(dva, dva)
+
+    pg = view(stack.vuser, k+1:k+func.ngen)
+    dpg = pg .- func.pgen_ref
+    q3 = 0.5 * func.πp * dot(dpg, dpg)
+
+    output[1] = q1 + q2 + q3
+    return
+end
+
+function adjoint!(func::TrackingCost, ∂stack, stack, ∂v)
+    k = func.k
+    vm = view(stack.vmag, 1:func.nbus)
+    ∂vm = view(∂stack.vmag, 1:func.nbus)
+    ∂vm .+= func.πv .* (vm .- func.vmag_ref) .* ∂v[1]
+
+    va = view(stack.vang, 1:func.nbus)
+    ∂va = view(∂stack.vang, 1:func.nbus)
+    ∂va .+= func.πa .* (va .- func.vang_ref) .* ∂v[1]
+
+    pg = view(stack.vuser, k+1:k+func.ngen)
+    shift = k
+    for i in 1:k
+        ∂pg = view(∂stack.vuser, shift+1:shift+func.ngen)
+        ∂pg .+= func.πp .* (pg .- func.pgen_ref) .* ∂v[1]
+        shift += func.ngen
+    end
+    return
+end
+
diff --git a/src/Polar/second_order.jl b/src/Polar/second_order.jl
index 286595b7..48629a2a 100644
--- a/src/Polar/second_order.jl
+++ b/src/Polar/second_order.jl
@@ -47,10 +47,10 @@ function HessianProd(polar::PolarForm{T, VI, VT, MT}, func::AutoDiff.AbstractExp
     t1s{N} = ForwardDiff.Dual{Nothing,Float64, N} where N
     VD = A{t1s{1}}
 
-    stack = NetworkStack(nbus, ngen, nlines, VT, VD)
+    stack = NetworkStack(nbus, ngen, nlines, polar.ncustoms, VT, VD)
     init!(polar, stack)
 
-    ∂stack = NetworkStack(nbus, ngen, nlines, VT, VD)
+    ∂stack = NetworkStack(nbus, ngen, nlines, polar.ncustoms, VT, VD)
 
     t1sF = zeros(Float64, n_cons) |> VD
     adj_t1sF = similar(t1sF)
@@ -93,6 +93,30 @@ function _hessian_sparsity(polar::AbstractPolarFormulation, func)
     return matpower_hessian(polar, func, V, y)
 end
 
+function _hessian_sparsity(polar::PolarFormRecourse, func)
+    m = length(func)
+    nbus, ngen = get(polar, PS.NumberOfBuses()), get(polar, PS.NumberOfGenerators())
+    gen = 2nbus+1:2nbus+ngen
+    Vre = Float64[i for i in 1:nbus]
+    Vim = Float64[i for i in nbus+1:2*nbus]
+    V = Vre .+ im .* Vim
+    y = rand(m)
+    H = matpower_hessian(polar, func, V, y)
+    n = size(H, 1)
+    @assert n == 2nbus + ngen
+    # Hessian is dimensioned without additional variables, pad it
+    H21 = spzeros(n, ngen+1)
+    H22 = [
+           ones(1, 1)   ones(1, ngen)
+           ones(ngen, 1) H[gen, gen]
+    ]
+
+    return [
+        H    H21
+        H21' H22
+    ]::SparseMatrixCSC
+end
+
 function _get_hessian_colors(polar::AbstractPolarFormulation, func::AutoDiff.AbstractExpression, map::Vector{Int})
     H = _hessian_sparsity(polar, func)::SparseMatrixCSC
     Hsub = H[map, map] # reorder
@@ -113,7 +137,7 @@ struct FullHessian{Model, Func, Stack, VD, SMT, VI} <: AutoDiff.AbstractFullHess
     H::SMT
 end
 
-function FullHessian(polar::PolarForm{T, VI, VT, MT}, func::AutoDiff.AbstractExpression, map::Vector{Int}) where {T, VI, VT, MT}
+function FullHessian(polar::AbstractPolarFormulation{T, VI, VT, MT}, func::AutoDiff.AbstractExpression, map::Vector{Int}) where {T, VI, VT, MT}
     (SMT, A) = get_jacobian_types(polar.device)
 
     pf = polar.network
@@ -133,10 +157,10 @@ function FullHessian(polar::PolarForm{T, VI, VT, MT}, func::AutoDiff.AbstractExp
     H = H_host |> SMT
 
     # Structures
-    stack = NetworkStack(nbus, ngen, nlines, VT, VD)
+    stack = NetworkStack(nbus, ngen, nlines, polar.ncustoms, VT, VD)
     init!(polar, stack)
 
-    ∂stack = NetworkStack(nbus, ngen, nlines, VT, VD)
+    ∂stack = NetworkStack(nbus, ngen, nlines, polar.ncustoms, VT, VD)
     t1sF = zeros(Float64, n_cons) |> VD
     adj_t1sF = similar(t1sF)
 
@@ -204,7 +228,7 @@ function hessian_arrowhead_sparsity(H, nx, nu, nscen)
 end
 
 function ArrowheadHessian(
-    polar::BlockPolarForm{T, VI, VT, MT},
+    polar::AbstractPolarFormulation{T, VI, VT, MT},
     func::AutoDiff.AbstractExpression,
     X::AbstractVariable,
 ) where {T, VI, VT, MT}
@@ -265,10 +289,10 @@ function ArrowheadHessian(
     H = sparse(i_coo, j_coo, ones(length(i_coo)), ntot, ntot) |> SMT
 
     # Structures
-    stack = NetworkStack(nbus, ngen, nlines, k, VT, VD)
+    stack = NetworkStack(nbus, ngen, nlines, polar.ncustoms, k, VT, VD)
     init!(polar, stack)
 
-    ∂stack = NetworkStack(nbus, ngen, nlines, k, VT, VD)
+    ∂stack = NetworkStack(nbus, ngen, nlines, polar.ncustoms, k, VT, VD)
     t1sF = zeros(Float64, n_cons) |> VD
     adj_t1sF = similar(t1sF)
 
diff --git a/src/Polar/stacks.jl b/src/Polar/stacks.jl
index 5c38a475..25cc799c 100644
--- a/src/Polar/stacks.jl
+++ b/src/Polar/stacks.jl
@@ -61,6 +61,7 @@ struct NetworkStack{VT,VD,NT} <: AbstractNetworkStack{VT}
     vmag::VD # voltage magnitudes
     vang::VD # voltage angles
     pgen::VD # active power generations
+    vuser::VD # custom variables defined by user
     # INTERMEDIATE
     ψ::VD    # nonlinear basis ψ(vmag, vang)
     intermediate::NT
@@ -71,8 +72,8 @@ struct NetworkStack{VT,VD,NT} <: AbstractNetworkStack{VT}
     nblocks::Int
 end
 
-function NetworkStack(nbus::Int, ngen::Int, nlines::Int, k::Int, VT::Type, VD::Type)
-    m = (2 * nbus + ngen) * k
+function NetworkStack(nbus::Int, ngen::Int, nlines::Int, nuser::Int, k::Int, VT::Type, VD::Type)
+    m = (2 * nbus + ngen) * k + nuser
     input = VD(undef, m) ; fill!(input, 0.0)
     # Wrap directly array x to avoid dealing with views
     p0 = pointer(input)
@@ -81,6 +82,12 @@ function NetworkStack(nbus::Int, ngen::Int, nlines::Int, k::Int, VT::Type, VD::T
     vang = unsafe_wrap(VD, p1, k * nbus)
     p2 = pointer(input, 2*k*nbus+1)
     pgen = unsafe_wrap(VD, p2, k * ngen)
+    p3 = pointer(input, 2*k*nbus+k*ngen+1)
+    vuser = if nuser > 0
+        unsafe_wrap(VD, p3, nuser)
+    else
+        similar(VD, 0)
+    end
 
     # Basis function
     ψ = VD(undef, k*(2*nlines + nbus)) ; fill!(ψ, 0.0)
@@ -107,17 +114,17 @@ function NetworkStack(nbus::Int, ngen::Int, nlines::Int, k::Int, VT::Type, VD::T
     p1 = pointer(params, k*nbus+1)
     qload = unsafe_wrap(VT, p1, k*nbus)
 
-    return NetworkStack(input, vmag, vang, pgen, ψ, intermediate, params, pload, qload, k)
+    return NetworkStack(input, vmag, vang, pgen, vuser, ψ, intermediate, params, pload, qload, k)
 end
-function NetworkStack(nbus::Int, ngen::Int, nlines::Int, VT::Type, VD::Type)
-    return NetworkStack(nbus, ngen, nlines, 1, VT, VD)
+function NetworkStack(nbus::Int, ngen::Int, nlines::Int, ncustoms::Int, VT::Type, VD::Type)
+    return NetworkStack(nbus, ngen, nlines, ncustoms, 1, VT, VD)
 end
 
 function NetworkStack(polar::AbstractPolarFormulation{T,VI,VT,MT}) where {T,VI,VT,MT}
     nbus = get(polar, PS.NumberOfBuses())
     ngen = get(polar, PS.NumberOfGenerators())
     nlines = get(polar, PS.NumberOfLines())
-    stack = NetworkStack(nbus, ngen, nlines, nblocks(polar), VT, VT)
+    stack = NetworkStack(nbus, ngen, nlines, polar.ncustoms, nblocks(polar), VT, VT)
     init!(polar, stack)
     return stack
 end
diff --git a/src/PowerSystem/PowerSystem.jl b/src/PowerSystem/PowerSystem.jl
index d94a2cd8..933a9b6b 100644
--- a/src/PowerSystem/PowerSystem.jl
+++ b/src/PowerSystem/PowerSystem.jl
@@ -259,6 +259,23 @@ struct Branches{T}
     to_buses::Vector{Int}
 end
 
+number_lines(b::Branches) = length(b.Yff)
+
+function drop_line(b::Branches, line_id::Int)
+    nl = number_lines(b)
+    @assert 1 <= line_id <= nl
+    # Screening
+    lines_on = ones(nl)
+    lines_on[line_id] = 0
+    return Branches{Complex{Float64}}(
+        b.Yff .* lines_on,
+        b.Yft .* lines_on,
+        b.Ytf .* lines_on,
+        b.Ytt .* lines_on,
+        b.from_buses, b.to_buses,
+    )
+end
+
 # Utils
 function get_bus_id_to_indexes(bus)
     BUS_I = IndexSet.idx_bus()[1]
diff --git a/src/PowerSystem/power_network.jl b/src/PowerSystem/power_network.jl
index e74721b3..a4ba77ef 100644
--- a/src/PowerSystem/power_network.jl
+++ b/src/PowerSystem/power_network.jl
@@ -119,6 +119,16 @@ function PowerNetwork(datafile::String; options...)
     return PowerNetwork(data; options...)
 end
 
+function PowerNetwork(network::PowerNetwork)
+    data = Dict{String, Array}()
+    data["bus"] = copy(network.buses)
+    data["branch"] = copy(network.branches)
+    data["gen"] = copy(network.generators)
+    data["baseMVA"] = Float64[network.baseMVA]
+    data["cost"] = copy(network.costs)
+    return PowerNetwork(data)
+end
+
 # Wrap ExaData and PGLIB's artifacts
 @enum(PowerNetworkLibrary,
     EXADATA=1,
@@ -350,11 +360,16 @@ function get_costs_coefficients(pf::PowerNetwork)
     return coefficients
 end
 
-function get_basis_matrix(pf::PowerNetwork)
-    nb = pf.nbus
-    nl = size(pf.branches, 1)
-    Yff, Yft, Ytf, Ytt = pf.lines.Yff, pf.lines.Yft, pf.lines.Ytf, pf.lines.Ytt
-    f, t = pf.lines.from_buses, pf.lines.to_buses
+function get_basis_matrix(pf::PowerNetwork; remove_line::Int=0)
+    nb, nl = pf.nbus, number_lines(pf.lines)
+    # Screening
+    lines = if 1 <= remove_line <= nl
+        drop_line(pf.lines, remove_line)
+    else
+        pf.lines
+    end
+
+    f, t = lines.from_buses, lines.to_buses
 
     Cf = sparse(f, 1:nl, ones(nl), nb, nl)       # connection matrix for line & from buses
     Ct = sparse(t, 1:nl, ones(nl), nb, nl)       # connection matrix for line & to buses
@@ -363,25 +378,29 @@ function get_basis_matrix(pf::PowerNetwork)
     Ysh = sparse(1:nb, 1:nb, ysh, nb, nb)
 
     # Build matrix
-    Yc = Cf * Diagonal(Yft) + Ct * Diagonal(Ytf)
-    Ys = Cf * Diagonal(Yft) - Ct * Diagonal(Ytf)
-    Yd = Cf * Diagonal(Yff) * Cf' + Ct * Diagonal(Ytt) * Ct' + Ysh
+    Yc = Cf * Diagonal(lines.Yft) + Ct * Diagonal(lines.Ytf)
+    Ys = Cf * Diagonal(lines.Yft) - Ct * Diagonal(lines.Ytf)
+    Yd = Cf * Diagonal(lines.Yff) * Cf' + Ct * Diagonal(lines.Ytt) * Ct' + Ysh
 
-    return [-real(Yc) -imag(Ys) -real(Yd);
+    return [-real(Yc)  -imag(Ys) -real(Yd);
              imag(Yc)  -real(Ys)  imag(Yd)]
 end
 
-function get_line_flow_matrices(pf::PowerNetwork)
-    nb = pf.nbus
-    nl = size(pf.branches, 1)
-    Yff, Yft, Ytf, Ytt = pf.lines.Yff, pf.lines.Yft, pf.lines.Ytf, pf.lines.Ytt
+function get_line_flow_matrices(pf::PowerNetwork; remove_line::Int=0)
+    nb, nl = pf.nbus, number_lines(pf.lines)
+    # Screening
+    lines = if 1 <= remove_line <= nl
+        drop_line(pf.lines, remove_line)
+    else
+        pf.lines
+    end
 
-    yff = Diagonal(Yff)
-    yft = Diagonal(Yft)
-    ytf = Diagonal(Ytf)
-    ytt = Diagonal(Ytt)
+    yff = Diagonal(lines.Yff)
+    yft = Diagonal(lines.Yft)
+    ytf = Diagonal(lines.Ytf)
+    ytt = Diagonal(lines.Ytt)
 
-    f, t = pf.lines.from_buses, pf.lines.to_buses
+    f, t = lines.from_buses, lines.to_buses
 
     Cf = sparse(f, 1:nl, ones(nl), nb, nl)       # connection matrix for line & from buses
     Ct = sparse(t, 1:nl, ones(nl), nb, nl)       # connection matrix for line & to buses
diff --git a/src/cuda_wrapper.jl b/src/cuda_wrapper.jl
index 8f1a0842..263b3460 100644
--- a/src/cuda_wrapper.jl
+++ b/src/cuda_wrapper.jl
@@ -4,11 +4,16 @@ import CUDA.CUBLAS
 import CUDA.CUSPARSE: CuSparseMatrixCSR, CuSparseMatrixCSC
 import CUDA.CUSOLVER
 
-function PolarForm(pf::PS.PowerNetwork, device::CUDABackend)
-    return PolarForm{Float64, CuVector{Int}, CuVector{Float64}, CuMatrix{Float64}}(pf, device)
+function PolarForm(pf::PS.PowerNetwork, device::CUDABackend, ncustoms::Int=0)
+    return PolarForm{Float64, CuVector{Int}, CuVector{Float64}, CuMatrix{Float64}}(pf, device, ncustoms)
 end
-function BlockPolarForm(pf::PS.PowerNetwork, device::CUDABackend, k::Int)
-    return BlockPolarForm{Float64, CuVector{Int}, CuVector{Float64}, CuMatrix{Float64}}(pf, device, k)
+function BlockPolarForm(pf::PS.PowerNetwork, device::CUDABackend, k::Int, ncustoms::Int=0)
+    return BlockPolarForm{Float64, CuVector{Int}, CuVector{Float64}, CuMatrix{Float64}}(pf, device, k, ncustoms)
+end
+function PolarFormRecourse(pf::PS.PowerNetwork, device::CUDABackend, k::Int)
+    ngen = PS.get(pf, PS.NumberOfGenerators())
+    ncustoms = (ngen + 1) * k
+    return PolarFormRecourse{Float64, CuVector{Int}, CuVector{Float64}, CuMatrix{Float64}}(pf, device, k, ncustoms)
 end
 
 default_sparse_matrix(::CUDABackend) = CuSparseMatrixCSR{Float64, Int32}
diff --git a/test/Polar/TestPolarForm.jl b/test/Polar/TestPolarForm.jl
index 58142a52..cae549fc 100644
--- a/test/Polar/TestPolarForm.jl
+++ b/test/Polar/TestPolarForm.jl
@@ -18,6 +18,7 @@ const LS = LinearSolvers
 include("api.jl")
 include("first_order.jl")
 include("second_order.jl")
+include("recourse.jl")
 
 function myisless(a, b)
     h_a = a |> Array
@@ -48,9 +49,6 @@ function runtests(case, device, AT)
         test_polar_stack(polar, device, AT)
         test_polar_constraints(polar, device, AT)
         test_polar_powerflow(polar, device, AT)
-        test_polar_blockstack(polar, device, AT)
-        test_polar_blk_expressions(polar, device, AT)
-        test_block_powerflow(polar, device, AT)
     end
 
     @testset "PolarForm AutoDiff (first-order)" begin
@@ -58,14 +56,34 @@ function runtests(case, device, AT)
         test_constraints_adjoint(polar, device, AT)
         test_full_space_jacobian(polar, device, AT)
         test_reduced_gradient(polar, device, AT)
-        test_block_jacobian(polar, device, AT)
     end
 
     @testset "PolarForm AutoDiff (second-order)" begin
         test_hessprod_with_finitediff(polar, device, AT)
         test_full_space_hessian(polar, device, AT)
+    end
+
+    @testset "BlockPolarForm" begin
+        test_block_stack(polar, device, AT)
+        test_block_expressions(polar, device, AT)
+        test_block_powerflow(polar, device, AT)
+        test_block_jacobian(polar, device, AT)
         test_block_hessian(polar, device, AT)
     end
+
+    @testset "Contingency" begin
+        test_contingency_powerflow(polar, device, AT)
+    end
+
+    @testset "PolarFormRecourse" begin
+        test_recourse_expression(polar, device, AT)
+        test_recourse_powerflow(polar, device, AT)
+        if isa(device, CPU)
+            test_recourse_jacobian(polar, device, AT)
+            test_recourse_hessian(polar, device, AT)
+            test_recourse_block_hessian(polar, device, AT)
+        end
+    end
 end
 
 end
diff --git a/test/Polar/api.jl b/test/Polar/api.jl
index 13756191..f693eea1 100644
--- a/test/Polar/api.jl
+++ b/test/Polar/api.jl
@@ -41,7 +41,7 @@ function test_polar_stack(polar, device, M)
     return nothing
 end
 
-function test_polar_blockstack(polar, device, M)
+function test_block_stack(polar, device, M)
     nblocks = 2
     blk_polar = ExaPF.BlockPolarForm(polar, nblocks)
     ngen = get(polar, PS.NumberOfGenerators())
@@ -212,7 +212,7 @@ function test_polar_powerflow(polar, device, M)
     end
 end
 
-function test_polar_blk_expressions(polar, device, M)
+function test_block_expressions(polar, device, M)
     nblocks = 2
     blk_polar = ExaPF.BlockPolarForm(polar, nblocks)
     stack = ExaPF.NetworkStack(polar)
@@ -289,3 +289,52 @@ function test_block_powerflow(polar, device, M)
     @test convergence.norm_residuals < pf_solver.tol
 end
 
+function test_contingency_powerflow(polar, device, M)
+    nbus = get(polar, PS.NumberOfBuses())
+    pf_solver = NewtonRaphson(tol=1e-10)
+
+    contingencies = ExaPF.LineContingency[
+        ExaPF.LineContingency(2),
+        ExaPF.LineContingency(8),
+    ]
+    nblocks = length(contingencies) + 1
+
+    blk_polar = ExaPF.BlockPolarForm(polar, nblocks)
+    blk_stack = ExaPF.NetworkStack(blk_polar)
+
+    pf = ExaPF.PowerFlowBalance(blk_polar, contingencies) ∘ ExaPF.PolarBasis(blk_polar)
+
+    blk_jac = ExaPF.ArrowheadJacobian(blk_polar, pf, State())
+    ExaPF.set_params!(blk_jac, blk_stack)
+    ExaPF.jacobian!(blk_jac, blk_stack)
+
+    convergence = ExaPF.nlsolve!(
+        pf_solver,
+        blk_jac,
+        blk_stack;
+    )
+    @test convergence.has_converged
+    @test convergence.norm_residuals < pf_solver.tol
+
+    # Compare with references
+    for (k, contingency) in enumerate(contingencies)
+        line = contingency.line_id
+        network = PS.PowerNetwork(polar.network)
+        # Drop line manually inside model's data
+        network.lines.Yff[line] = 0.0im
+        network.lines.Yft[line] = 0.0im
+        network.lines.Ytf[line] = 0.0im
+        network.lines.Ytt[line] = 0.0im
+        #
+        post_model = ExaPF.PolarForm(network, device)
+
+        stack = ExaPF.NetworkStack(post_model)
+        conv = ExaPF.run_pf(post_model, stack)
+        @test conv.has_converged
+
+        vref = Array(stack.vmag)
+        vmag = Array(blk_stack.vmag[k * nbus + 1:(k+1)*nbus])
+        @test vref ≈ vmag
+    end
+end
+
diff --git a/test/Polar/recourse.jl b/test/Polar/recourse.jl
new file mode 100644
index 00000000..7ae098c4
--- /dev/null
+++ b/test/Polar/recourse.jl
@@ -0,0 +1,196 @@
+
+function test_recourse_powerflow(polar, device, M)
+    k = 2
+    polar_ext = ExaPF.PolarFormRecourse(polar, k)
+    stack = ExaPF.NetworkStack(polar_ext)
+
+    pf_recourse = ExaPF.PowerFlowRecourse(polar_ext) ∘ ExaPF.PolarBasis(polar_ext)
+    jac_recourse = ExaPF.ArrowheadJacobian(polar_ext, pf_recourse, State())
+    ExaPF.set_params!(jac_recourse, stack)
+    ExaPF.jacobian!(jac_recourse, stack)
+
+    pf_solver = NewtonRaphson()
+    convergence = ExaPF.nlsolve!(
+        pf_solver,
+        jac_recourse,
+        stack,
+    )
+    @test convergence.has_converged
+    @test convergence.norm_residuals < pf_solver.tol
+
+    residual = pf_recourse(stack)
+    @test norm(residual) < pf_solver.tol
+    return
+end
+
+function test_recourse_expression(polar, device, M)
+    k = 2
+    polar_ext = ExaPF.PolarFormRecourse(polar, k)
+    @test ExaPF.nblocks(polar_ext) == k
+
+    stack = ExaPF.NetworkStack(polar_ext)
+
+    ngen = PS.get(polar, PS.NumberOfGenerators())
+    @test polar_ext.ncustoms == k * (ngen + 1)
+    @test length(stack.vuser) == k * (ngen + 1)
+
+    for expr in [
+        ExaPF.PowerFlowRecourse,
+        ExaPF.ReactivePowerBounds,
+    ]
+        ev = expr(polar_ext) ∘ ExaPF.PolarBasis(polar_ext)
+        res = ev(stack)
+        @test isa(res, M)
+        @test length(res) == length(ev)
+
+        g_min, g_max = ExaPF.bounds(polar_ext, ev)
+        @test length(g_min) == length(ev)
+        @test length(g_max) == length(ev)
+        @test isa(g_min, M)
+        @test isa(g_max, M)
+        @test myisless(g_min, g_max)
+    end
+    return
+end
+
+function test_recourse_jacobian(polar, device, M)
+    k = 2
+    polar_ext = ExaPF.PolarFormRecourse(polar, k)
+    stack = ExaPF.NetworkStack(polar_ext)
+    ∂stack = ExaPF.NetworkStack(polar_ext)
+    stack_fd = ExaPF.NetworkStack(polar_ext)
+
+    nu = ExaPF.number(polar_ext, Control())
+    mapx = ExaPF.mapping(polar_ext, State(), k)
+    mapu = ExaPF.mapping(polar_ext, Control(), k)
+    mapxu = ExaPF.mapping(polar_ext, AllVariables(), k)
+
+    for expr in [
+        ExaPF.PowerFlowRecourse,
+        ExaPF.ReactivePowerBounds,
+    ]
+        ev = expr(polar_ext) ∘ ExaPF.PolarBasis(polar_ext)
+        m = length(ev)
+
+        # Compute ref with finite-diff
+        function _fd_func(x)
+            stack_fd.input .= x
+            return ev(stack_fd)
+        end
+        x0 = copy(stack.input)
+        Jd = FiniteDiff.finite_difference_jacobian(_fd_func, x0) |> Array
+        Jd_x = Jd[:, mapx]
+        Jd_u = sum(Jd[:, mapu[1+(i-1)*nu:i*nu]] for i in 1:k)
+        Jd_xu = [Jd_x Jd_u]
+
+        # / State
+        jac_x = ExaPF.ArrowheadJacobian(polar_ext, ev, State())
+        Jx = ExaPF.jacobian!(jac_x, stack) |> SparseMatrixCSC
+        @test Jx ≈ Jd_x rtol=1e-5
+
+        # / Control
+        jac_u = ExaPF.ArrowheadJacobian(polar_ext, ev, Control())
+        Ju = ExaPF.jacobian!(jac_u, stack) |> SparseMatrixCSC
+        @test Ju ≈ Jd_u rtol=1e-5
+
+        # / all
+        jac_xu = ExaPF.ArrowheadJacobian(polar_ext, ev, AllVariables())
+        Jxu = ExaPF.jacobian!(jac_xu, stack) |> SparseMatrixCSC
+        @test Jxu ≈ Jd_xu rtol=1e-5
+
+        # Test adjoint
+        tgt_h = rand(m)
+        tgt = tgt_h |> M
+        empty!(∂stack)
+        ExaPF.adjoint!(ev, ∂stack, stack, tgt)
+        # Compare with finite diff
+        function _fd_adj(x)
+            stack_fd.input .= x
+            return dot(tgt_h, ev(stack_fd))
+        end
+        adj_fd = FiniteDiff.finite_difference_gradient(_fd_adj, x0)
+        @test myisapprox(∂stack.input[mapxu], adj_fd[mapxu], rtol=1e-6)
+        @test myisapprox(∂stack.input[mapxu], Jd[:, mapxu]' * tgt_h, rtol=1e-6)
+    end
+    return
+end
+
+function test_recourse_hessian(polar, device, M)
+    k = 1
+    polar_ext = ExaPF.PolarFormRecourse(polar, k)
+    stack = ExaPF.NetworkStack(polar_ext)
+    stack_fd = ExaPF.NetworkStack(polar_ext)
+    ∂stack = ExaPF.NetworkStack(polar_ext)
+
+    basis  = ExaPF.PolarBasis(polar_ext)
+    mapxu = ExaPF.mapping(polar_ext, AllVariables())
+
+    constraints = [
+        ExaPF.QuadraticCost(polar_ext),
+        ExaPF.PowerFlowRecourse(polar_ext),
+        ExaPF.ReactivePowerBounds(polar_ext),
+        ExaPF.LineFlows(polar_ext),
+    ]
+    ev = ExaPF.MultiExpressions(constraints) ∘ basis
+
+    m = length(ev)
+    c = zeros(m)
+    y = rand(m)
+
+    # Evaluate Hessian
+    hess = ExaPF.FullHessian(polar_ext, ev, mapxu)
+    H = ExaPF.hessian!(hess, stack, y)
+
+    # Evaluate Hessian with finite-diff
+    function _fd_hess(x)
+        stack_fd.input[mapxu] .= x
+        return dot(y, ev(stack_fd))
+    end
+    x0 = stack.input[mapxu]
+    Hd = FiniteDiff.finite_difference_hessian(_fd_hess, x0)
+    @test myisapprox(H, Hd, rtol=1e-5)
+    return
+end
+
+function test_recourse_block_hessian(polar, device, M)
+    k = 2
+    polar_ext = ExaPF.PolarFormRecourse(polar, k)
+    stack = ExaPF.NetworkStack(polar_ext)
+    stack_fd = ExaPF.NetworkStack(polar_ext)
+    ∂stack = ExaPF.NetworkStack(polar_ext)
+
+    basis  = ExaPF.PolarBasis(polar_ext)
+    mapxu = ExaPF.mapping(polar_ext, State(), k)
+
+    constraints = [
+        ExaPF.PowerFlowRecourse(polar_ext),
+        ExaPF.ReactivePowerBounds(polar_ext),
+        ExaPF.LineFlows(polar_ext),
+    ]
+    ev = ExaPF.MultiExpressions(constraints) ∘ basis
+
+    m = length(ev)
+    c = zeros(m)
+    y = rand(m)
+
+    # Evaluate Hessian
+    hess = ExaPF.ArrowheadHessian(polar_ext, ev, AllVariables())
+    H = ExaPF.hessian!(hess, stack, y)
+
+    # Evaluate Hessian with finite-diff
+    function _fd_grad(x)
+        stack_fd.input[mapxu] .= x
+        ev(c, stack_fd)
+        empty!(∂stack)
+        ExaPF.adjoint!(ev, ∂stack, stack_fd, y)
+        return ∂stack.input[mapxu]
+    end
+    x0 = stack.input[mapxu]
+    Hd = FiniteDiff.finite_difference_jacobian(_fd_grad, x0)
+
+    # Test evaluation matches with finite-diff
+    # TODO
+    # @test myisapprox(H, Hd, rtol=1e-5)
+    return H, Hd
+end
+