Skip to content

simplify hessian matrix structure initialization#349

Merged
jd-lara merged 6 commits into
mainfrom
lk/homotopy-hessian-jtj
May 14, 2026
Merged

simplify hessian matrix structure initialization#349
jd-lara merged 6 commits into
mainfrom
lk/homotopy-hessian-jtj

Conversation

@luke-kiernan
Copy link
Copy Markdown
Collaborator

@luke-kiernan luke-kiernan commented Apr 7, 2026

Currently there's some graph stuff for figuring out the sparse structure of the hessian of the homotopy function. I did some profiling a while back and found that simply computing J^T J is faster. Both give the same results: we want (2-or-fewer)-hop neighbors. Also using J^T J's structure generalizes better: e.g. proportional-to-headroom slack.

Copy link
Copy Markdown
Contributor

Copilot AI left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Pull request overview

This PR simplifies initialization of the robust-homotopy Hessian sparse structure by deriving it directly from the Jacobian product J' * J (instead of building the structure via a graph-based neighbor traversal), motivated by profiling/performance.

Changes:

  • Remove the custom graph-based _create_hessian_matrix_structure implementation.
  • Allocate the Hessian sparsity pattern by temporarily filling J.Jv.nzval, computing Hv = J.Jv' * J.Jv, and zeroing stored values afterward.

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

Comment thread src/RobustHomotopy/homotopy_hessian.jl Outdated
Comment thread src/RobustHomotopy/homotopy_hessian.jl Outdated
@github-actions
Copy link
Copy Markdown
Contributor

github-actions Bot commented Apr 7, 2026

Performance Results

Precompile Time

Main This Branch Delta
2.486 s 2.368 s -4.7%

Solve Time

Test Main This Branch Delta
matpower_ACTIVSg10k_sys-NewtonRaphsonACPowerFlow First Solve 12.012 s 11.872 s -1.2%
matpower_ACTIVSg10k_sys-NewtonRaphsonACPowerFlow Second Solve 81.7 ms 538.6 ms +559.4%
matpower_ACTIVSg10k_sys-RobustHomotopyPowerFlow First Solve 16.114 s 15.749 s -2.3%
matpower_ACTIVSg10k_sys-RobustHomotopyPowerFlow Second Solve 9.623 s 9.261 s -3.8%
matpower_ACTIVSg10k_sys-NewtonRaphsonACPowerFlow(iwamoto) First Solve 393.3 ms 202.4 ms -48.5%
matpower_ACTIVSg10k_sys-NewtonRaphsonACPowerFlow(iwamoto) Second Solve 80.2 ms 87.0 ms +8.5%
matpower_ACTIVSg10k_sys-TrustRegionACPowerFlow(iwamoto) First Solve 2.158 s 2.441 s +13.1%
matpower_ACTIVSg10k_sys-TrustRegionACPowerFlow(iwamoto) Second Solve 463.7 ms 85.5 ms -81.6%
matpower_ACTIVSg10k_sys-DCPowerFlow First Solve 5.127 s 5.088 s -0.8%
matpower_ACTIVSg10k_sys-DCPowerFlow Second Solve 14.6 ms 14.7 ms +0.8%
matpower_ACTIVSg10k_sys-PTDFDCPowerFlow First Solve 1.719 s 1.709 s -0.6%
matpower_ACTIVSg10k_sys-PTDFDCPowerFlow Second Solve 65.4 ms 65.4 ms +0.1%
matpower_ACTIVSg10k_sys-vPTDFDCPowerFlow First Solve 5.943 s 5.925 s -0.3%
matpower_ACTIVSg10k_sys-vPTDFDCPowerFlow Second Solve 3.561 s 3.542 s -0.5%

luke-kiernan and others added 2 commits April 7, 2026 14:51
@luke-kiernan
Copy link
Copy Markdown
Collaborator Author

luke-kiernan commented Apr 8, 2026

Copilot was right — verified that ac_power_flow_jacobian.jl:331-332 sets the LCC angle-constraint diagonals to 1.0 at structure creation, and _update_jacobian_matrix_values! never rewrites them, so the old nonzeros(J.Jv) .= 0.0 would corrupt them on any system with LCC. Fixed via save/restore of J.Jv's nzval in d9b58e6.

edit: turns out the RH method as written isn't compatible with LCCs anyway...but it's a one-time cost, so I'll leave the save-and-restore.

@luke-kiernan
Copy link
Copy Markdown
Collaborator Author

luke-kiernan commented Apr 8, 2026

Turns out these safeguards aren't strictly necessary: at the moment, H = J' * J keeps structural zeros. However, nothing in the SparseArrays docs dictates that behavior, so best not to rely on it.

LCCs add state variables to the Jacobian that the homotopy Hessian
formulation does not account for, producing a dimension mismatch
deep inside the solve. Fail fast at HomotopyHessian construction
with a clear error message.

Co-Authored-By: Claude Opus 4.7 (1M context) <[email protected]>
@luke-kiernan
Copy link
Copy Markdown
Collaborator Author

I added a "no LCCs" guard on the RH method path and a basic test that exercises said guard.

@luke-kiernan luke-kiernan requested a review from jd-lara May 14, 2026 20:29
@jd-lara
Copy link
Copy Markdown
Member

jd-lara commented May 14, 2026

I added a "no LCCs" guard on the RH method path and a basic test that exercises said guard.

The problem is going to the that in real systems like WECC and the EI there will be LCCs in the data.

@luke-kiernan
Copy link
Copy Markdown
Collaborator Author

luke-kiernan commented May 14, 2026

The problem is going to the that in real systems like WECC and the EI there will be LCCs in the data.

Sure, I can work on an LCC-compatible robust homotopy method. But I'd still vote for merging this: throwing a human-readable error is improvement over "dimension mismatch."

edit: there's a couple other things in this PR too, besides the LCC error

@jd-lara
Copy link
Copy Markdown
Member

jd-lara commented May 14, 2026

The problem is going to the that in real systems like WECC and the EI there will be LCCs in the data.

Sure, I can work on an LCC-compatible robust homotopy method. But I'd still vote for merging this: throwing a human-readable error is improvement over "dimension mismatch."

edit: there's a couple other things in this PR too, besides the LCC error

Yeah, I don't think is blocking. Just nothing that it will be a consideration

@jd-lara jd-lara merged commit dbbf578 into main May 14, 2026
7 checks passed
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

3 participants