Batch power flow
ExaPF provides a way to evaluate the expressions by blocks, opening the way to introduce more parallelism in the code.
BlockPolarForm
We recall that a given NetworkStack stack stores the different variables and parameters (power generations, voltages, loads) required to evaluate the power flow model.
stack = ExaPF.NetworkStack(polar);21-elements NetworkStack{Vector{Float64}}The variables are stored in the field stack.input, the parameters in the field stack.params. The parameters encode the active pd and reactive loads qd at all buses in the network, such that
nbus = ExaPF.get(polar, PS.NumberOfBuses());
pd = stack.params[1:nbus]
qd = stack.params[nbus+1:2*nbus]9-element Vector{Float64}:
0.0
0.0
0.0
0.0
0.3
0.0
0.35
0.0
0.5By default, a NetworkStack stores one set of loads $p_0$.
Suppose now we want to evaluate the model associated with the polar formulation for $N$ different set of parameters (=scenarios) $p_1, \cdots, p_N$. ExaPF allows to streamline the polar formulation with a BlockPolarForm structure:
nscen = 10;
blk_polar = ExaPF.BlockPolarForm(polar, nscen)10-BlockPolar formulation (instantiated on device CPU(false))
Network characteristics:
#buses: 9 (#slack: 1 #PV: 2 #PQ: 6)
#generators: 3
#lines: 9
giving a mathematical formulation with:
#controls: 5
#states : 14Then, ExaPF can also instantiate a NetworkStack object, with the memory required to store the variables of the different scenarios:
blk_stack = ExaPF.NetworkStack(blk_polar)210-elements NetworkStack{Vector{Float64}}We can pass the scenarios manually using the function set_params!:
ploads = rand(nbus, nscen);
qloads = rand(nbus, nscen);
ExaPF.set_params!(blk_stack, ploads, qloads)90-element Vector{Float64}:
0.030498754705436948
0.7972832375644416
0.42036867854961235
0.2445133620603247
0.6496158927308912
0.5116658390581503
0.404459007639985
0.899372753729086
0.21903480541958043
0.6165836415810982
⋮
0.21279943172895588
0.12259419188181753
0.8262082589393523
0.6738699674853635
0.16453130808943162
0.1480083789204356
0.28332357811047426
0.3748667165011048
0.8960630860500634The structure blk_stack stores $N$ different realizations for the variables stored in the field input (vmag, vang and pgen). By default, the initial values are set according to the values specified in blk_polar (usually defined when importing the data from the instance file):
reshape(blk_stack.vmag, nbus, nscen)9×10 Matrix{Float64}:
1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0
1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0
1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0
1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0
1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0
1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0
1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0
1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0
1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0Only the parameters are varying according to the scenarios we passed as input in the constructor:
reshape(blk_stack.pload, nbus, nscen)9×10 Matrix{Float64}:
0.156243 0.841015 0.830368 0.280914 … 0.0975969 0.915708 0.567575
0.596326 0.362657 0.887451 0.608459 0.896085 0.0866123 0.783787
0.725615 0.0319489 0.64645 0.677092 0.880479 0.168619 0.713891
0.897547 0.848453 0.351371 0.498256 0.0735903 0.694895 0.59484
0.732423 0.973632 0.689133 0.206782 0.357461 0.697483 0.365112
0.84924 0.745735 0.949638 0.770137 … 0.500139 0.713763 0.984281
0.944556 0.370409 0.398403 0.861146 0.791993 0.88803 0.0611511
0.744242 0.0598381 0.639591 0.635662 0.640238 0.947311 0.372265
0.0878342 0.577902 0.446943 0.157641 0.164826 0.911853 0.336649Evaluate expressions in block
ExaPF takes advantage of the block structure when using a BlockPolarForm.
As an example, suppose we want to evaluate the power flow balances in block form with a PowerFlowBalance expression:
powerflow = ExaPF.PowerFlowBalance(blk_polar) ∘ ExaPF.PolarBasis(blk_polar);ExaPF.ComposedExpressions{ExaPF.PolarBasis{Vector{Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}}, ExaPF.PowerFlowBalance{Vector{Float64}, SparseArrays.SparseMatrixCSC{Float64, Int64}}}(PolarBasis (AbstractExpression), PowerFlowBalance (AbstractExpression))A block evaluation takes as input the NetworkStack blk_stack structure:
m = div(length(powerflow), nscen);
blk_output = zeros(m * nscen);
powerflow(blk_output, blk_stack); # inplace evaluation
reshape(blk_output, m, nscen)14×10 Matrix{Float64}:
-1.03367 -1.26734 -0.742549 … -0.733915 -1.54339 -0.846213
-0.124385 -0.818051 -0.20355 0.0304794 -0.681381 -0.136109
0.897547 0.848453 0.351371 0.0735903 0.694895 0.59484
0.732423 0.973632 0.689133 0.357461 0.697483 0.365112
0.84924 0.745735 0.949638 0.500139 0.713763 0.984281
0.944556 0.370409 0.398403 … 0.791993 0.88803 0.0611511
0.744242 0.0598381 0.639591 0.640238 0.947311 0.372265
0.0878342 0.577902 0.446943 0.164826 0.911853 0.336649
0.0775134 -0.13659 -0.12913 -7.12952e-5 -0.140896 0.50687
0.391616 0.392255 0.200366 0.0179658 0.398802 -0.0934687
0.228166 0.177695 0.301178 … -0.0605341 0.564021 -0.135492
0.225459 0.476406 0.714421 -0.107935 0.667575 0.104324
0.671873 -0.0546481 0.108339 0.582329 0.25019 0.147367
-0.0219652 0.446682 -0.190313 0.685207 -0.130531 0.655063We get $N$ different results for the power flow balance equations, depending on which scenario we are on.
Solve power flow in block on the CPU
Once the different structures used for block evaluation instantiated, one is able to solve the power flow in block on the CPU using the same function nlsolve!. The block Jacobian is evaluated with automatic differentiation using a ArrowheadJacobian structure:
blk_jx = ExaPF.ArrowheadJacobian(blk_polar, powerflow, State());
blk_jx.J140×140 SparseArrays.SparseMatrixCSC{Float64, Int64} with 820 stored entries:
⎡⡱⣮⡲⣞⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎤
⎢⣸⢮⣻⣾⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⡱⣮⡲⣞⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⣸⢮⣻⣾⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⡵⣯⡶⣝⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⣜⢯⡻⣮⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡵⣯⡶⣝⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣜⢯⡻⣮⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡵⣯⡶⣝⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣜⢯⡻⣮⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡵⣯⣺⣜⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣚⢾⡻⣮⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡵⣯⣺⣜⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣚⢾⡻⣮⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡵⣯⣺⣜⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣚⢾⡻⣮⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡵⣯⡺⣕⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢞⢮⡿⣯⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡵⣯⡺⣕⎥
⎣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢞⢮⡿⣯⎦We notice that the ArrowheadJacobian computes the resulting Jacobian as a block diagonal matrix. The ArrowheadJacobian has a slightly different behavior than its classical counterpart AutoDiff.Jacobian, in the sense that one has to pass the parameters manually to initiate internally the dual numbers:
ExaPF.set_params!(blk_jx, blk_stack);
ExaPF.jacobian!(blk_jx, blk_stack);140×140 SparseArrays.SparseMatrixCSC{Float64, Int64} with 820 stored entries:
⎡⡱⣮⡲⣞⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎤
⎢⣸⢮⣻⣾⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⡱⣮⡲⣞⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⣸⢮⣻⣾⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⡵⣯⡶⣝⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⣜⢯⡻⣮⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡵⣯⡶⣝⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣜⢯⡻⣮⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡵⣯⡶⣝⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣜⢯⡻⣮⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡵⣯⣺⣜⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣚⢾⡻⣮⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡵⣯⣺⣜⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣚⢾⡻⣮⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡵⣯⣺⣜⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣚⢾⡻⣮⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡵⣯⡺⣕⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢞⢮⡿⣯⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡵⣯⡺⣕⎥
⎣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢞⢮⡿⣯⎦As soon as the blk_jx initialized, we can solve the power flow equations in block as
conv = ExaPF.nlsolve!(
NewtonRaphson(verbose=2),
blk_jx,
blk_stack;
)Power flow has converged: true
* #iterations: 4
* Time Jacobian (s) ........: 0.0004
* Time linear solver (s) ...: 0.0001
* update (s) ............: 0.0001
* ldiv (s) ..............: 0.0000
* Time total (s) ...........: 0.0009
At the solution, we get different values for the voltage magnitudes at the PQ nodes:
reshape(blk_stack.vmag, nbus, nscen)9×10 Matrix{Float64}:
1.0 1.0 1.0 1.0 … 1.0 1.0 1.0
1.0 1.0 1.0 1.0 1.0 1.0 1.0
1.0 1.0 1.0 1.0 1.0 1.0 1.0
0.951858 0.959992 0.982733 0.938902 0.96507 0.968818 0.950918
0.901319 0.917828 0.942606 0.93944 0.957697 0.917554 0.960086
0.942211 0.961329 0.94828 0.967046 … 0.975752 0.9349 0.981508
0.9186 0.940148 0.917821 0.940863 0.958752 0.90356 0.967187
0.933614 0.968099 0.959164 0.933488 0.948122 0.946842 0.967972
0.943275 0.92586 0.980698 0.896571 0.912839 0.95912 0.911801Solve power flow in batch on the GPU
When the BlockPolarForm model is instantiated on the GPU, the expressions are evaluated in batch. The syntax to solve the power flow equations is exactly the same as on the CPU, using cusolverRF to solve the different linear systems:
using CUDA, CUSOLVERRF
polar_gpu = ExaPF.load_polar("case9.m", CUDABackend());
blk_polar_gpu = ExaPF.BlockPolarForm(polar_gpu, nscen); # load model on GPU
blk_stack_gpu = ExaPF.NetworkStack(blk_polar_gpu);
ExaPF.set_params!(blk_stack_gpu, ploads, qloads);
powerflow_gpu = ExaPF.PowerFlowBalance(blk_polar_gpu) ∘ ExaPF.PolarBasis(blk_polar_gpu);
blk_jx_gpu = ExaPF.ArrowheadJacobian(blk_polar_gpu, powerflow_gpu, State());
ExaPF.set_params!(blk_jx_gpu, blk_stack_gpu);
ExaPF.jacobian!(blk_jx_gpu, blk_stack_gpu);
rf_fac = CUSOLVERRF.RFLU(blk_jx_gpu.J)
rf_solver = LS.DirectSolver(rf_fac)
conv = ExaPF.nlsolve!(
NewtonRaphson(verbose=2),
blk_jx_gpu,
blk_stack_gpu;
linear_solver=rf_solver,
)Power flow has converged: true
* #iterations: 4
* Time Jacobian (s) ........: 0.0023
* Time linear solver (s) ...: 0.1157
* update (s) ............: 0.0585
* ldiv (s) ..............: 0.0572
* Time total (s) ...........: 0.1219