Evaluator: compute reduced Hessian on GPU #118
Conversation
* new function to compute Hessian of objective with AutoDiff, via Hessian-vector product * rename `cost_production` as `objective` * rename `\partial cost()` as `adjoint_objective!`
Allow to evaluate Hessian on the GPU directly ReducedSpaceEvaluator: * implement hessprod! with AutoDiff * implement hessprod_lagrangian_penalty! with AutoDiff * replace MATPOWER's Jacobian with AutoDiff's Jacobian * remove old attribute `autodiff` in ReducedSpaceEvaluator * improve type inference for AutoDiff * remove method `get(nlp, AutoDiffBackend())`, which was ambiguous * remove Hessian and Jacobian caches in ReducedSpaceEvaluator Polar: * remove init_autodiff_factory in Polar/ * add a new ConstantJacobian structure for linear function * add a new function `is_linear` to state whether a constraint is linear * add new storage types `JacobianStorage` and `HessianStorage`
| @@ -55,9 +55,9 @@ powerflow_solver = NewtonRaphson(tol=ntol) | |||
| nlp = ExaPF.ReducedSpaceEvaluator(polar, x0, u0; | |||
| linear_solver=algo, powerflow_solver=powerflow_solver) | |||
| convergence = ExaPF.update!(nlp, u0) | |||
| ExaPF.reset!(nlp) | |||
There was a problem hiding this comment.
Lately I tried to get better at this. I have no idea why I spontaneously spread this around the code. Must be VisualCode's fault :-P.
There was a problem hiding this comment.
No doubt that this is VisualCode's fault!
| end | ||
|
|
||
| function ReducedSpaceEvaluator( | ||
| model, x, u; | ||
| constraints=Function[voltage_magnitude_constraints, active_power_constraints, reactive_power_constraints], | ||
| linear_solver=DirectSolver(), | ||
| powerflow_solver=NewtonRaphson(tol=1e-12), | ||
| want_jacobian=true, | ||
| want_hessian=true, |
| reduced_gradient!(nlp, jv, jvx, jvx, jvu) | ||
| return | ||
| end | ||
|
|
||
| function ojtprod!(nlp::ReducedSpaceEvaluator, jv, u, σ, v) |
There was a problem hiding this comment.
Is this the adjoint of the objective? Could be called gradient.
There was a problem hiding this comment.
This is not exactly the adjoint of the objective. It's the adjoint of the vector-value function x -> [f(x); h(x)] with f the objective, h the nonlinear constraints
| # ψ = -(∇gₓ' \ (∇²fx .+ ∇²gx)) | ||
| if isa(hessvec, CUDA.CuArray) | ||
| ∇gₓ = nlp.state_jacobian.Jx.J | ||
| ∇gT = LinearSolvers.get_transpose(nlp.linear_solver, ∇gₓ) |
There was a problem hiding this comment.
Aha, I should have looked at that.
There was a problem hiding this comment.
I am thinking of a cleaner implementation for all our linear operators. Let's discuss that today :-)
| # Jacobian | ||
| ∇gₓ = nlp.autodiff.Jgₓ.J::SpMT | ||
| ∇gᵤ = nlp.autodiff.Jgᵤ.J::SpMT | ||
| # Two vector products |
There was a problem hiding this comment.
This is like a cache for preallocation?
| @@ -125,6 +125,36 @@ function AutoDiff.jacobian!(polar::PolarForm, jac::AutoDiff.Jacobian, buffer) | |||
| return jac.J | |||
| end | |||
|
|
|||
| # Handle properly constant Jacobian case | |||
| function AutoDiff.ConstantJacobian(polar::PolarForm, func::Function, variable::Union{State,Control}) | |||
There was a problem hiding this comment.
Nice. Have you run into this case?
There was a problem hiding this comment.
Indeed! The constraints w.r.t. voltage magnitude are all linear, leading to a constant Jacobian.
* change signature of `hessian_prod_objective` * remove deprecated functions * rename adjoint_objective! as gradient_objective!
Allow to evaluate Hessian on the GPU directly
ReducedSpaceEvaluator:
autodiffin ReducedSpaceEvaluatorget(nlp, AutoDiffBackend()), whichwas ambiguous
Polar:
via Hessian-vector product
is_linearto state whether a constraint is linearJacobianStorageandHessianStoragecost_productionasobjective\partial cost()asadjoint_objective!