ModelingToolkit.jl, an IR and Compiler for Scientific Models
juliadiffeq.org
juliadiffeq.org
As someone who spends more time on symbolic mathematics than numerical mathematics, I’ve often wished that Julia’s symbolic libraries had the same sort of energy behind them that the numerical libraries had. But it’s great to see that a lot of the work put into the numerical libraries is relevant to the symbolic ones.
(this is a compliment! it's great stuff)
When you model the transmission grid, you can solve for the flows on the equipment (lines and transformers) where the inputs are certain equipment parameters as well as generation and load values. This is called powerflow or loadflow and the traditional Newton-Raphson method represents the equations in Jacobian form where the sensitivities represent how real and reactive power is impacted by the bus voltages and bus phase angles. This solution takes voltage into account (pretty important) and forms the basis of many computationally intensive operations such as contingency analysis (also known as N-1). You basically determine before the analysis which equipment you want to monitor and then you have a list of contingencies that you want to guard against. The grid itself is a pretty massive NxN matrix that is thankfully both sparse and symmetric. When you simulate a contingency, you take it out of service in the model and then re-solve the AC powerflow. For each of those solves, you check that the bus voltages are ok and that the equipment is within the rated limits. Violations are reported to the user. There are other calcs that are also important and that use the Jacobian. I think one reason this area hasn't received much GPU attention is that the problem is easy enough to parallelize and most vendors do so now, so a server with a decent amount of cores can whittle away at even large grids with thousands of contingencies in a reasonable amount of time. Less important studies that need to be quick and robust use the DC powerflow which is linearized and doesn't show voltage impact.
There are many entities that monitor this in real-time, but transmission planners work with text file based models and will run millions of contingencies with computer clusters in a batch mode and dump the results in Excel or a SQL database. I'm not sure if a high end GPU would help there or if it is just a throw more CPU cores kind of problem.
How much would it be worth? It is tough to say. I'd have to know how well GPUs perform here. I'm guessing it is not as good, but I'm not knowledgeable here. If you had a product that bundled GPU, GPU-based power systems analysis software and a client that would run on desktop with an API, AND it performed significantly faster than the competition, you might have some customers. That is a ton of upfront work though.
Does anyone (I assume the poster I originally replied to) know which situations are even worth it to is solve the Jacobian on a GPU?
One of the questions I have about the Julia diffeq module is how many simultaneous equations can it take? Is it thousands?
If you question is, how far can ModelingToolkit be pushed? That's unknown. We made sure it's quite an efficient representation, but we haven't stress tested it to see the memory consumption at a million operations yet. Since it's able to generate thousands of operations in symbolic Jacobian inversions, it easy goes to thousands and tens of thousands of equations. How much larger? Who knows. What we can say is that its efficiency is more similar to SymEngine's symbolic types, so if SymEngine handles the millions of equations range (on a computer with enough RAM, symbolic representations may use up quite a bit here!) than we should be able to.
I'm looking forward to starting to play with these new libraries. Is there a more focused place to talk about Julia applied to numerical work (linear algebra and nonlinear equation solvers).
Drop me a line if you want: Certhas@gmail.com
I'll send you an email sometime in the near future for some discussion. If you're using Julia for grid work, you're probably either one of the firms listed on their success stories, or one of the national labs. Either of which sounds interesting from one engineer to another.
But what can we do with ModelingToolkit to start pushing more into this domain? Well, here's a few things:
1) We plan to have new compilation targets, kind of like the SArray one that was shown, which target parallelism. For example, for large Jacobians, we want a target so you can say "generate the code that fills the Jacobian optimally with 16 threads", and then what it will do is split up the Jacobian into 16 parts with approximately equally many operations and thread between those. Since it's all scalar operations then this should be something that gets as efficient as you'd want to write by hand! Then a similar compilation target can be made for the GPU, where it would need to know the number of GPU threads you'd be using and it would split apart that Jacobian filling and write a gigantic `if threadid() == ...` setup, and use CUDAnative.jl to compile this to a .ptx kernel so that all of those threads are parallelized. Both of these trivially extend to sparse matrices doing what I showed in the writing. So filling sparse matrices is something we can and plan to have nice compilation targets for.
2) We can also have compilation targets for `J*v` operations, so that way it can accelerate Newton-Krylov solvers, which is probably the better way of using sparse matrices on the GPU since it would avoid sparse factorizations.
3) We can also implement a generic sparse factorization in Julia, to then build analytical solutions for sparse factorizations, run Rewrite.jl simplification passes on these, and generate the parallelized kernels for them. This would then be something that does sparse factorizations on the GPU with parallelized code but does not need to use CUSOLVER, and so that would be quite ideal. The issue is that it would require writing generic sparse factorization tools, since Julia currently relies on SuiteSparse for this with no generic fallback. This is not an easy task, but I think Yingbo will be up for it in a few years since we have a plan to not use standard sparse linear algebra routines in Julia's differential equation suites because there is a way to speed up operations specifically in this domain by using a new factorization, and so if we go this direction our first task will be to learn how the current sparse LU is implemented.
I'm really excited for what y'all have coming in the future and will stay tuned. There is some literature on Newton-Krylov solvers for my domain (power systems) so I'll have to brush up on that math.
Help me understand one thing. I suppose if we could one day compile Julia to a binary and sell software (just curious...I don't sell software today), would you have to secure a license for suitesparse? Is that another benefit from rewriting?