Is your code public?
Yeah that’s not enough for sparse diff. (That’s barely sparse ). Maybe the decompression could be optimized a bit more though.
What do you consider sparse? the number of structural non-zeros is roughly 3% for the IEEE118 case:
(NonZeros, Entries) = (1506, 55696)
And the percentage only goes down for larger cases. For the 2383 case we have:
(NonZeros, Entries) = (30636, 22714756) (MaxColor, Size) = (20, 4766)
That’s roughly 0.15% of structural non-zeros.
SymbolicAD seems to leverage this structure properly though.
I just uploaded it to Github. Just run
test.jl and that’s it.
EDIT: I forgot to mention that right now I’m only computing the PF Hessians and Jacobians. I will compute those terms for the constraints as well, but I haven’t implemented that yet.
Perhaps one other clarification. We’re only computing the jacobian and hessian of nonlinear (i.e., not linear or quadratic) constraints and objectives. We don’t compute it for linear or quadratic terms since they get handled by the solver.
Does that include power flow constraints? they are general nonlinear in polar coordinates.
Yes, what @odow is indicating is that these are the only equations that
SymbolicAD is managing for building Jacobian and Hessian. The other parts of the AC-OPF problem are quadratic and have a different interface to the solver due to their simplified Jac/Hess structure.
It has been a couple of months since I made this post and I am planning re-run these test with latest versions of all packages. Please let me know of any specific package improvements that you think could impact this analysis.
The new Optimization.jl sparsity handling is probably worth a try.
Additionally, if @YingboMa gets around to it then it might be interesting to try tearing
I literally had a draft open with reply to your original post that we added sparsity support in Optimization.jl so you might want to give it a try
But we haven’t added documentation for it yet, the way to do it is pass
AutoModelingToolkit(true, true) (the two fields are boolean for enabling sparsity for objective and constraint respectively).
I finished implementing sparsity and symbolic differentiation support in Nonconvex. Your rosetta-opf example works locally. I will make a release and update the docs shortly.
@mohamed82008, I tried with Nonconvex v1.0.3 but I am still getting this issue of zygote crashing, could you elaborate on how I should change these calls to use your symbolic diff approach?
model = Nonconvex.DictModel() ... x0 = NonconvexCore.getinit(model) r = Nonconvex.optimize(model, IpoptAlg(), x0)
First you will need Nonconvex v1.0.4. Then you can do:
first_order = true options = IpoptOptions(; first_order) sp_model = sparsify(model, hessian = !first_order) r = Nonconvex.optimize(sp_model, IpoptAlg(), x0; options)
sparsify will use Symbolics.jl to get the sparsity pattern of the gradients/Jacobians/Hessians and then use
ForwardDiff for the differentiation. Alternatively, you can use:
first_order = true options = IpoptOptions(; first_order) sym_model = symbolify(model, hessian = !first_order, sparse = true, simplify = true) r = Nonconvex.optimize(sym_model, IpoptAlg(), x0; options)
which will use a
Symbolics-only approach where the gradients/Jacobians/Hessians are obtained using Symbolics.jl. The
sparse kwarg is used to tell Symbolics to return a sparse matrix/vector for the jacobian/Hessian/gradient.
Note that these are new features so please report any issues in the repo.
As promised, below is a table of updated results using the latest versions of all packages. Note that GalacticOptim has been renamed Optimization.
The substantive changes follow the discussion above. In particular, NonConvex and Optimization have been updated with new AD systems based on Symbolics and now support sparse Jacobian and Hessian computations. For NonConvex, this allowed me to finish building the NonConvex OPF model and confirm that the model is correct on small networks. For Optimization, the new sparse AD system based on ModelingTookit enabled it to solve models with 73 network buses (800+ variables, 900+ constraints), which is a 2x scalability increase from what was possible with the previous dense ForwardDiff approach.
For general awareness,
- This Zygote issue has thus far blocked using Zygote for AD on these models
- These upstream issues  appear to be the root cause blocking scaling of the ModelingTookKit approach to larger problems. It looks like it will be resolved in future Julia versions (testing was done on Julia v1.7.2).
 ADNLPModels v0.3.3
[b6b21f68] Ipopt v1.0.2
[4076af6c] JuMP v1.1.1
[961ee093] ModelingToolkit v8.14.0
[f4238b75] NLPModelsIpopt v0.9.1
[01bcebdf] Nonconvex v2.0.3
[bf347577] NonconvexIpopt v0.4.1
[429524aa] Optim v1.7.0
[7f7a1694] Optimization v3.6.1
[fd9f6733] OptimizationMOI v0.1.1
[c36e90e8] PowerModels v0.19.5
Ipopt was run with the default linear solver, MUMPS 5.4.1.
Are you including the compilation time here? Yeah the RGF issues are one thing, and the other is we need to release SymbolicIR.
MTK also doesn’t do
structural_simplify on OptimizationProblem right now, which can be a major win because then it solves a much smaller problem. @YingboMa was working on that.
Good question. In these runtime numbers I try to avoid capturing artifacts from one-time overheads, such as JIT time. In my experience, with sufficient effort, these can be avoided in production deployment settings.
In all of these implementations before I start collecting data
solve_opf("$(@__DIR__)/data/pglib_opf_case5_pjm.m") is executed first, to compile the relevant code execution paths. Then
solve_opf(...) is run again on the file where runtime data is being collected.
You can have a look at exactly how the data collection is done in this script. It is not immediately obvious that pre-compilation is done there, but looking into the specific implementation files you will find it, for example, rosetta-opf/optimization.jl at cli · lanl-ansi/rosetta-opf · GitHub
Cool, looking forward to seeing the next iteration.
Thanks for this excellent source of test cases for Nonconvex. For Nonconvex, I would be interested in the results using ReverseDiff as well as SparseDiffTools instead of Symbolics.
For SparseDiffTools, just drop the
symbolic = true:
result = Nonconvex.optimize( model, IpoptAlg(), NonconvexCore.getinit(model); options = IpoptOptions(; first_order=false, sparse=true), )
For ReverseDiff, you can do:
using AbstractDifferentiation, ReverseDiff rd_model = abstractdiffy(model, AD.ReverseDiffBackend()) result = Nonconvex.optimize( rd_model, IpoptAlg(), NonconvexCore.getinit(model); options = IpoptOptions(; first_order=false), )
I don’t make use of tape compilation yet in Nonconvex but ReverseDiff is often quite fast already without compiling a tape for each function.
Also how can I run all these problems for Nonconvex locally?