The ExaData artifact contains test cases relevant to the Exascale Computing Project. It is built from the git repository available at ExaData. Apart from the standard MATPOWER files it additionally contains demand scenarios and contingencies used in multiperiod security constrained optimal power flow settings.
Settings
This document was generated with Documenter.jl version 1.0.1 on Friday 22 September 2023. Using Julia version 1.9.3.
+ ${display_result} +
+ExaPF.jl is a package to solve the power flow problem on upcoming exascale architectures. The code has been designed to be:
ExaPF implements a vectorized modeler for power systems, which allows to manipulate basic expressions. All expressions are fully differentiable : their first and second-order derivatives can be extracted efficiently using automatic differentiation. In addition, we leverage the packages CUDA.jl and KernelAbstractions.jl to make ExaPF portable across GPU architectures.
This research was supported by the Exascale Computing Project (17-SC-20-SC), a joint project of the U.S. Department of Energy’s Office of Science and National Nuclear Security Administration, responsible for delivering a capable exascale ecosystem, including software, applications, and hardware technology, to support the nation’s exascale computing imperative.
Settings
This document was generated with Documenter.jl version 1.0.1 on Friday 22 September 2023. Using Julia version 1.9.3.
ExaPF.AutoDiff.AbstractStack — TypeAbstractStack{VT}Abstract variable storing the inputs and the intermediate values in the expression tree.
ExaPF.AutoDiff.AbstractExpression — TypeAbstractExpressionAbstract type for differentiable function $f(x)$. Any AbstractExpression implements two functions: a forward mode to evaluate $f(x)$, and an adjoint to evaluate $∂f(x)$.
Forward mode
The direct evaluation of the function $f$ is implemented as
(expr::AbstractExpression)(output::VT, stack::AbstractStack{VT}) where VT<:AbstractArray
+the input being specified in stack, the results being stored in the array output.
Reverse mode
The adjoint of the function is specified by the function adjoint!, with the signature:
adjoint!(expr::AbstractExpression, ∂stack::AbstractStack{VT}, stack::AbstractStack{VT}, ̄v::VT) where VT<:AbstractArray
+The variable stack stores the result of the direct evaluation, and is not modified in adjoint!. The results are stored inside the adjoint stack ∂stack.
ExaPF.AutoDiff.adjoint! — Functionadjoint!(expr::AbstractExpression, ∂stack::AbstractStack{VT}, stack::AbstractStack{VT}, ̄v::VT) where VT<:AbstractArrayCompute the adjoint of the AbstractExpression expr with relation to the adjoint vector ̄v. The results are stored in the adjoint stack ∂stack. The variable stack stores the result of a previous direct evaluation, and is not modified in adjoint!.
ExaPF.AutoDiff.AbstractJacobian — TypeAbstractJacobianAutomatic differentiation for the compressed Jacobian of any nonlinear constraint $h(x)$.
ExaPF.AutoDiff.AbstractHessianProd — TypeAbstractHessianProdReturns the adjoint-Hessian-vector product $λ^⊤ H v$ of any nonlinear constraint $h(x)$.
ExaPF.AutoDiff.AbstractFullHessian — TypeAbstractHessianProdFull sparse Hessian $H$ of any nonlinear constraint $h(x)$.
ExaPF.AutoDiff.seed! — Functionseed!(
+ H::AbstractHessianProd,
+ v::AbstractVector{T},
+) where {T}Seed the duals with v to compute the Hessian vector product $λ^⊤ H v$.
ExaPF.AutoDiff.seed_coloring! — Functionseed_coloring!(
+ M::Union{AbstractJacobian, AbstractFullHessian}
+ coloring::AbstractVector,
+)Seed the duals with the coloring based seeds to compute the Jacobian or Hessian $M$.
ExaPF.AutoDiff.partials! — Functionpartials!(jac::AbstractJacobian)Extract partials from Jacobian jac in jac.J.
partials!(hess::AbstractFullHessian)Extract partials from Hessian hess into hess.H.
ExaPF.AutoDiff.set_value! — Functionset_value!(
+ jac,
+ primals::AbstractVector{T}
+) where {T}Set values of ForwardDiff.Dual numbers in jac to primals.
Settings
This document was generated with Documenter.jl version 1.0.1 on Friday 22 September 2023. Using Julia version 1.9.3.
ExaPF.AbstractVariable — TypeAbstractVariableVariables corresponding to a particular formulation.
ExaPF.AbstractFormulation — TypeAbstractFormulationInterface between the data and the mathemical formulation.
ExaPF.State — TypeState <: AbstractVariableAll variables $x$ depending on the variables Control $u$ through the non-linear equation $g(x, u) = 0$.
ExaPF.Control — TypeControl <: AbstractVariableIndependent variables $u$ used in the reduced-space formulation.
ExaPF.PolarForm — TypePolarForm{T, IT, VT, MT} <: AbstractPolarFormulationImplement the polar formulation associated to the network's equations.
Wrap a PS.PowerNetwork network to load the data on the target device (CPU() and CUDABackend() are currently supported).
Example
julia> const PS = ExaPF.PowerSystem;
+
+julia> network_data = PS.load_case("case9.m");
+
+julia> polar = PolarForm(network_data, ExaPF.CPU())
+Polar 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 : 14
+ExaPF.BlockPolarForm — TypeBlockPolarForm{T, IT, VT, MT} <: AbstractFormulationBlock polar formulation: duplicates k different polar models to evaluate them in parallel.
ExaPF.load_polar — Functionload_polar(case, device=CPU(); dir=PS.EXADATA)Load a PolarForm instance from the specified benchmark library dir on the target device (default is CPU). ExaPF uses two different benchmark libraries: MATPOWER (dir=EXADATA) and PGLIB-OPF (dir=PGLIB).
Examples
julia> polar = ExaPF.load_polar("case9")
+Polar 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 : 14
+ExaPF.NetworkStack — TypeNetworkStack{VT,VD,MT} <: AbstractNetworkStack{VT}
+NetworkStack(polar::PolarForm)
+NetworkStack(nbus::Int, ngen::Int, nlines::Int, VT::Type)Store the variables associated to the polar formulation. The variables are stored in the field input, ordered as follows
input = [vmag ; vang ; pgen]The object stores also intermediate variables needed in the expression tree, such as the LKMR basis ψ.
Notes
The NetworkStack can be instantiated on the host or on the target device.
Examples
julia> polar = ExaPF.load_polar("case9");
+
+julia> stack = ExaPF.NetworkStack(polar)
+21-elements NetworkStack{Vector{Float64}}
+
+julia> stack.vmag
+9-element Vector{Float64}:
+ 1.0
+ 1.0
+ 1.0
+ 1.0
+ 1.0
+ 1.0
+ 1.0
+ 1.0
+ 1.0
+ExaPF.init! — Functioninit!(polar::PolarForm, stack::NetworkStack)Set stack.input with the initial values specified in the base PS.PowerNetwork object.
The state and the control are defined as mapping:
ExaPF.mapping — Functionmapping(polar::PolarForm, ::Control)Return the mapping associated to the Control() in NetworkStack according to the polar formulation PolarForm.
Examples
julia> polar = ExaPF.load_polar("case9");
+
+julia> mapu = ExaPF.mapping(polar, Control())
+5-element Vector{Int64}:
+ 1
+ 2
+ 3
+ 20
+ 21
+mapping(polar::PolarForm, ::State)Return the mapping associated to the State() in NetworkStack according to the polar formulation PolarForm.
Examples
julia> polar = ExaPF.load_polar("case9");
+
+julia> mapu = ExaPF.mapping(polar, State())
+14-element Vector{Int64}:
+ 11
+ 12
+ 13
+ 14
+ 15
+ 16
+ 17
+ 18
+ 4
+ 5
+ 6
+ 7
+ 8
+ 9
+ExaPF.run_pf — Functionrun_pf(
+ polar::PolarForm, stack::NetworkStack;
+ rtol=1e-8, max_iter=20, verbose=0,
+)Solve the power flow equations $g(x, u) = 0$ w.r.t. the stack $x$, using the (NewtonRaphson algorithm. The initial state $x$ is specified implicitly inside stack, with the mapping mapping associated to the polar formulation. The object stack is modified inplace in the function.
The algorithm stops when a tolerance rtol or a maximum number of iterations maxiter is reached.
Arguments
polar::AbstractFormulation: formulation of the power flow equationstack::NetworkStack: initial values in the networkExamples
julia> polar = ExaPF.load_polar("case9");
+
+julia> stack = ExaPF.NetworkStack(polar);
+
+julia> conv = run_pf(polar, stack);
+
+julia> conv.has_converged
+true
+
+julia> conv.n_iterations
+4ExaPF.nlsolve! — Functionnlsolve!(
+ algo::NewtonRaphson,
+ jac::Jacobian,
+ stack::NetworkStack;
+ linear_solver=DirectSolver(jac.J),
+ nl_buffer=NLBuffer(size(jac, 2)),
+)Solve the nonlinear system of equations $g(x) = 0$ with a NewtonRaphson algorithm. At each iteration, we update the variable $x$ as
\[ x_{k+1} = x_{k} - (∇g_k)^{-1} g(x_k) +\]
till $\| g(x_k) \| < ε_{tol}$
In the implementation,
jac.func,stack::NetworkStack (with mapping jac.map),jac, with automatic differentiation.Note that stack is modified inplace during the iterations of algorithm.
The Jacobian jac should be instantied before calling this function. By default, the linear system $(∇g_k)^{-1} g(x_k)$ is solved using a LU factorization. You can specify a different linear solver by changing the optional argument linear_solver.
Arguments
algo::NewtonRaphon: Newton-Raphson object, storing the options of the algorithmjac::Jacobian: Stores the function $g$ and its Jacobian $∇g$. The Jacobian is updated with automatic differentiation.stack::NetworkStack: initial valueslinear_solver::AbstractLinearSolver: linear solver used to compute the Newton stepnl_buffer::NLBuffer: buffer storing the residual vector and the descent direction Δx. Can be reused to avoid unecessary allocations.Examples
julia> polar = ExaPF.load_polar("case9");
+
+julia> powerflow = ExaPF.PowerFlowBalance(polar) ∘ ExaPF.PolarBasis(polar);
+
+julia> jx = ExaPF.Jacobian(polar, powerflow, State());
+
+julia> stack = ExaPF.NetworkStack(polar);
+
+julia> conv = ExaPF.nlsolve!(NewtonRaphson(), jx, stack);
+
+julia> conv.has_converged
+true
+
+julia> conv.n_iterations
+4ExaPF.NewtonRaphson — TypeNewtonRaphson <: AbstractNonLinearSolverNewton-Raphson algorithm.
Attributes
maxiter::Int (default 20): maximum number of iterationstol::Float64 (default 1e-8): tolerance of the algorithmverbose::Int (default 0): verbosity levelThe different parts of the polar formulation are implemented in the following AbstractExpression:
ExaPF.PolarBasis — TypePolarBasis{VI, MT} <: AbstractExpression
+PolarBasis(polar::AbstractPolarFormulation)Implement the LKMR nonlinear basis. Takes as input the voltage magnitudes vmag and the voltage angles vang and returns
\[ \begin{aligned} + & \psi_\ell^C(v, \theta) = v^f v^t \cos(\theta_f - \theta_t) \quad \forall \ell = 1, \cdots, n_\ell \\ + & \psi_\ell^S(v, \theta) = v^f v^t \sin(\theta_f - \theta_t) \quad \forall \ell = 1, \cdots, n_\ell \\ + & \psi_k(v, \theta) = v_k^2 \quad \forall k = 1, \cdots, n_b + \end{aligned}\]
Dimension: 2 * n_lines + n_bus
Complexity
3 n_lines + n_bus mul, n_lines cos and n_lines sin
Examples
julia> polar = ExaPF.load_polar("case9");
+
+julia> stack = ExaPF.NetworkStack(polar);
+
+julia> basis = ExaPF.PolarBasis(polar)
+PolarBasis (AbstractExpression)
+
+julia> basis(stack)
+27-element Vector{Float64}:
+ 1.0
+ 1.0
+ 1.0
+ 1.0
+ 1.0
+ 1.0
+ 1.0
+ 1.0
+ 1.0
+ 0.0
+ ⋮
+ 1.0
+ 1.0
+ 1.0
+ 1.0
+ 1.0
+ 1.0
+ 1.0
+ 1.0
+ 1.0
+ExaPF.PowerFlowBalance — TypePowerFlowBalance{VT, MT}
+PowerFlowBalance(polar)Implement a subset of the power injection corresponding to $(p_{inj}^{pv}, p_{inj}^{pq}, q_{inj}^{pq})$. The function encodes the active balance equations at PV and PQ nodes, and the reactive balance equations at PQ nodes:
\[\begin{aligned} + p_i &= v_i \sum_{j}^{n} v_j (g_{ij}\cos{(\theta_i - \theta_j)} + b_{ij}\sin{(\theta_i - \theta_j})) \,, & + ∀ i ∈ \{PV, PQ\} \\ + q_i &= v_i \sum_{j}^{n} v_j (g_{ij}\sin{(\theta_i - \theta_j)} - b_{ij}\cos{(\theta_i - \theta_j})) \,. & + ∀ i ∈ \{PQ\} +\end{aligned}\]
Require composition with PolarBasis.
Dimension: n_pv + 2 * n_pq
Complexity
2 SpMV
Examples
julia> polar = ExaPF.load_polar("case9");
+
+julia> stack = ExaPF.NetworkStack(polar);
+
+julia> powerflow = ExaPF.PowerFlowBalance(polar) ∘ ExaPF.PolarBasis(polar);
+
+julia> round.(powerflow(stack); digits=6)
+14-element Vector{Float64}:
+ -1.63
+ -0.85
+ 0.0
+ 0.9
+ 0.0
+ 1.0
+ 0.0
+ 1.25
+ -0.167
+ 0.042
+ -0.2835
+ 0.171
+ -0.2275
+ 0.259
+
+julia> run_pf(polar, stack); # solve powerflow equations
+
+julia> isapprox(powerflow(stack), zeros(14); atol=1e-8)
+true
+ExaPF.VoltageMagnitudeBounds — TypeVoltageMagnitudeBoundsImplement the bounds on voltage magnitudes not taken into account in the bound constraints. In the reduced space, this is associated to the the voltage magnitudes at PQ nodes:
\[v_{pq}^♭ ≤ v_{pq} ≤ v_{pq}^♯ .\]
Dimension: n_pq
Complexity
1 copyto
Note
In the reduced space, the constraints on the voltage magnitudes at PV nodes $v_{pv}$ are taken into account when bounding the control $u$.
Examples
julia> polar = ExaPF.load_polar("case9");
+
+julia> stack = ExaPF.NetworkStack(polar);
+
+julia> voltage_pq = ExaPF.VoltageMagnitudeBounds(polar);
+
+julia> voltage_pq(stack)
+6-element Vector{Float64}:
+ 1.0
+ 1.0
+ 1.0
+ 1.0
+ 1.0
+ 1.0
+ExaPF.PowerGenerationBounds — TypePowerGenerationBounds{VT, MT}
+PowerGenerationBounds(polar)Constraints on the active power productions and on the reactive power productions that are not already taken into account in the bound constraints. In the reduced space, that amounts to
\[p_{g,ref}^♭ ≤ p_{g,ref} ≤ p_{g,ref}^♯ ; +C_g q_g^♭ ≤ C_g q_g ≤ C_g q_g^♯ .\]
Require composition with PolarBasis.
Dimension: n_pv + 2 n_ref
Complexity
1 copyto, 1 SpMV
Examples
julia> polar = ExaPF.load_polar("case9");
+
+julia> stack = ExaPF.NetworkStack(polar);
+
+julia> run_pf(polar, stack); # solve powerflow equations
+
+julia> power_generators = ExaPF.PowerGenerationBounds(polar) ∘ ExaPF.PolarBasis(polar);
+
+julia> round.(power_generators(stack); digits=6)
+4-element Vector{Float64}:
+ 0.719547
+ 0.24069
+ 0.144601
+ -0.03649
+ExaPF.LineFlows — TypeLineFlows{VT, MT}
+LineFlows(polar)Implement thermal limit constraints on the lines of the network.
Require composition with PolarBasis.
Dimension: 2 * n_lines
Complexity
4 SpMV, 4 * n_lines quadratic, 2 * n_lines add
Examples
julia> polar = ExaPF.load_polar("case9");
+
+julia> stack = ExaPF.NetworkStack(polar);
+
+julia> run_pf(polar, stack); # solve powerflow equations
+
+julia> line_flows = ExaPF.LineFlows(polar) ∘ ExaPF.PolarBasis(polar);
+
+julia> round.(line_flows(stack); digits=6)
+18-element Vector{Float64}:
+ 0.575679
+ 0.094457
+ 0.379983
+ 0.723832
+ 0.060169
+ 0.588673
+ 2.657418
+ 0.748943
+ 0.295351
+ 0.560817
+ 0.112095
+ 0.38625
+ 0.728726
+ 0.117191
+ 0.585164
+ 2.67781
+ 0.726668
+ 0.215497
+The production costs is given in the AbstractExpression CostFunction:
ExaPF.CostFunction — TypeCostFunction{VT, MT} <: AutoDiff.AbstractExpression
+CostFunction(polar)Implement the quadratic cost function for OPF
\[ ∑_{g=1}^{n_g} c_{2,g} p_g^2 + c_{1,g} p_g + c_{0,g}\]
Require composition with PolarBasis to evaluate the cost of the reference generator.
Dimension: 1
Complexity
1 SpMV, 1 sum
Examples
julia> polar = ExaPF.load_polar("case9");
+
+julia> stack = ExaPF.NetworkStack(polar);
+
+julia> cost = ExaPF.CostFunction(polar) ∘ ExaPF.PolarBasis(polar);
+
+julia> cost(stack)
+1-element Vector{Float64}:
+ 4509.0275
+The different expressions can be combined together in several different ways.
ExaPF.MultiExpressions — TypeMultiExpressions <: AbstractExpressionImplement expressions concatenation. Takes as input a vector of expressions [expr1,...,exprN] and concatenate them in a single expression mexpr, such that
mexpr(x) = [expr1(x) ; expr2(x) ; ... ; exprN(x)]
+ExaPF.ComposedExpressions — TypeComposedExpressions{Expr1<:PolarBasis, Expr2} <: AbstractExpressionImplement expression composition. Takes as input two expressions expr1 and expr2 and returns a composed expression cexpr such that ``` cexpr(x) = expr2 ∘ expr1(x)
Notes
Currently, only PolarBasis is supported for expr1.
Settings
This document was generated with Documenter.jl version 1.0.1 on Friday 22 September 2023. Using Julia version 1.9.3.
ExaPF allows to solve linear systems with either direct and indirect linear algebra, both on CPU and on GPU. To solve a linear system $Ax = b$, ExaPF uses the function ldiv!.
LinearAlgebra.ldiv! — Functionldiv!(solver, x, A, y)
+ldiv!(solver, x, y)solver::AbstractLinearSolver: linear solver to solve the systemx::AbstractVector: SolutionA::AbstractMatrix: Input matrixy::AbstractVector: RHSSolve the linear system $A x = y$ using the algorithm specified in solver. If A is not specified, the function will used directly the factorization stored inside solver.
ExaPF wraps UMFPACK (shipped with SuiteSparse.jl) on the CPU, and CUSPARSE on CUDA device.
ExaPF.LinearSolvers.DirectSolver — TypeDirectSolver <: AbstractLinearSolverSolve linear system $A x = y$ with direct linear algebra.
DirectSolver uses UMFPACK to solve the linear systemDirectSolver redirects the resolution to the method CUSOLVER.csrlsvqrExaPF.LinearSolvers.KrylovBICGSTAB — TypeKrylovBICGSTAB <: AbstractIterativeLinearSolver
+KrylovBICGSTAB(precond; verbose=0, rtol=1e-10, atol=1e-10)Wrap Krylov.jl's BICGSTAB algorithm to solve iteratively the linear system $A x = y$.
ExaPF.LinearSolvers.DQGMRES — TypeDQGMRES <: AbstractIterativeLinearSolver
+DQGMRES(precond; verbose=false, memory=4)Wrap Krylov.jl's DQGMRES algorithm to solve iteratively the linear system $A x = y$.
ExaPF.LinearSolvers.BICGSTAB — TypeBICGSTAB <: AbstractIterativeLinearSolver
+BICGSTAB(precond; maxiter=2_000, tol=1e-8, verbose=false)Custom BICGSTAB implementation to solve iteratively the linear system $A x = y$.
ExaPF.LinearSolvers.EigenBICGSTAB — TypeEigenBICGSTAB <: AbstractIterativeLinearSolver
+EigenBICGSTAB(precond; maxiter=2_000, tol=1e-8, verbose=false)Julia's port of Eigen's BICGSTAB to solve iteratively the linear system $A x = y$.
ExaPF.jl is shipped with a custom BICGSTAB implementation. However, we highly recommend to use KrylovBICGSTAB instead, which has proved to be more robust.
ExaPF.LinearSolvers.bicgstab — Functionbicgstab(A, b, P, xi;
+ tol=1e-8,
+ maxiter=size(A, 1),
+ verbose=false,
+ maxtol=1e20)BiCGSTAB implementation according to
Van der Vorst, Henk A. "Bi-CGSTAB: A fast and smoothly converging variant of Bi-CG for the solution of nonsymmetric linear systems." SIAM Journal on scientific and Statistical Computing 13, no. 2 (1992): 631-644.
Available linear solvers can be queried with
ExaPF.LinearSolvers.list_solvers — Functionlist_solvers(::KernelAbstractions.Device)List linear solvers available on current device.
To solve linear systems with iterative methods, ExaPF provides an implementation of a block-Jacobi preconditioner, portable on GPU.
ExaPF.LinearSolvers.AbstractPreconditioner — TypeAbstractPreconditionerPreconditioners for the iterative solvers mostly focused on GPUs
ExaPF.LinearSolvers.BlockJacobiPreconditioner — TypeBlockJacobiPreconditionerOverlapping-Schwarz preconditioner.
Attributes
nblocks::Int64: Number of partitions or blocks.blocksize::Int64: Size of each block.partitions::Vector{Vector{Int64}}:npart` partitions stored as listscupartitions: partitions transfered to the GPUlpartitions::Vector{Int64}`: Length of each partitions.culpartitions::Vector{Int64}`: Length of each partitions, on the GPU.blocks: Dense blocks of the block-Jacobicublocks: Js transfered to the GPUmap: The partitions as a mapping to construct viewscumap: cumap transferred to the GPU`part: Partitioning as output by Metiscupart: part transferred to the GPUExaPF.LinearSolvers.update — Functionfunction update(J::CuSparseMatrixCSR, p)Update the preconditioner p from the sparse Jacobian J in CSR format for the GPU
cuJs are filled from the sparse Jacobian Jp.P from the dense blocks cuJsfunction update(J::SparseMatrixCSC, p)Update the preconditioner p from the sparse Jacobian J in CSC format for the CPU
Note that this implements the same algorithm as for the GPU and becomes very slow on CPU with growing number of blocks.
Settings
This document was generated with Documenter.jl version 1.0.1 on Friday 22 September 2023. Using Julia version 1.9.3.
ExaPF.PowerSystem.AbstractPowerSystem — TypeAbstractPowerSystemFirst layer of the package. Store the topology of a given transmission network, including:
Data are imported either from a matpower file, or a PSSE file.
ExaPF.PowerSystem.PowerNetwork — TypePowerNetwork <: AbstractPowerSystemThis structure contains constant parameters that define the topology and physics of the power network.
The object PowerNetwork uses its own contiguous indexing for the buses. The indexing is independent from those specified in the Matpower or the PSSE input file. However, a correspondence between the two indexing (Input indexing to PowerNetwork indexing) is stored inside the attribute bus_to_indexes.
Note
The object PowerNetwork is created in the host memory. Use a AbstractFormulation object to move data to the target device.
ExaPF.PowerSystem.load_case — Functionload_case(case_name, lib::PowerNetworkLibrary=EXADATA)Convenient function to load a PowerNetwork instance from one of the benchmark library (dir=EXADATA for MATPOWER, dir=PGLIB for PGLIB-OPF). Default library is lib=EXADATA.
Examples
julia> net = PS.load_case("case118") # default is MATPOWER
+PowerNetwork object with:
+ Buses: 118 (Slack: 1. PV: 53. PQ: 64)
+ Generators: 54.
+
+julia> net = PS.load_case("case1354_pegase", PS.PGLIB)
+PowerNetwork object with:
+ Buses: 1354 (Slack: 1. PV: 259. PQ: 1094)
+ Generators: 260.ExaPF.PowerSystem.AbstractNetworkElement — TypeAbstractNetworkElementAbstraction for all physical elements being parts of a AbstractPowerSystem. Elements are divided in
Lines)Buses)Generators)List of elements:
ExaPF.PowerSystem.Buses — TypeBuses <: AbstractNetworkElementBuses of a transmission network.
ExaPF.PowerSystem.Lines — TypeLines <: AbstractNetworkElementLines of a transmission network.
ExaPF.PowerSystem.Generators — TypeGenerators <: AbstractElementGenerators in a transmission network
ExaPF.PowerSystem.AbstractNetworkAttribute — TypeAbstractNetworkAttributeAttribute of a AbstractPowerSystem.
List of attributes:
ExaPF.PowerSystem.NumberOfBuses — TypeNumberOfBuses <: AbstractNetworkAttributeNumber of buses in a AbstractPowerSystem.
ExaPF.PowerSystem.NumberOfLines — TypeNumberOfLines <: AbstractNetworkAttributeNumber of lines in a AbstractPowerSystem.
ExaPF.PowerSystem.NumberOfGenerators — TypeNumberOfGenerators <: AbstractNetworkAttributeNumber of generators in a AbstractPowerSystem.
ExaPF.PowerSystem.NumberOfPVBuses — TypeNumberOfPVBuses <: AbstractNetworkAttributeNumber of PV buses in a AbstractPowerSystem.
ExaPF.PowerSystem.NumberOfPQBuses — TypeNumberOfPQBuses <: AbstractNetworkAttributeNumber of PQ buses in a AbstractPowerSystem.
ExaPF.PowerSystem.NumberOfSlackBuses — TypeNumberOfSlackBuses <: AbstractNetworkAttributeNumber of slack buses in a AbstractPowerSystem.
ExaPF.PowerSystem.BaseMVA — TypeBaseMVA <: AbstractNetworkAttributeBase MVA of the network.
ExaPF.PowerSystem.BusAdmittanceMatrix — TypeBusAdmittanceMatrix <: AbstractNetworkAttributeBus admittance matrix associated with the topology of the network.
Query the indexing of the different elements in a given network:
ExaPF.PowerSystem.PVIndexes — TypePVIndexes <: AbstractIndexingIndexes of the PV buses in a AbstractPowerSystem.
ExaPF.PowerSystem.PQIndexes — TypePQIndexes <: AbstractIndexingIndexes of the PQ buses in a AbstractPowerSystem.
ExaPF.PowerSystem.SlackIndexes — TypeSlackIndexes <: AbstractIndexingIndexes of the slack buses in a AbstractPowerSystem.
ExaPF.PowerSystem.GeneratorIndexes — TypeGeneratorIndexes <: AbstractIndexingIndexes of the generators in a AbstractPowerSystem.
ExaPF.PowerSystem.AbstractNetworkValues — TypeAbstractNetworkValuesNumerical values attached to the different attributes of the network.
List of values:
ExaPF.PowerSystem.VoltageMagnitude — TypeVoltageMagnitude <: AbstractNetworkValuesMagnitude |v| of the voltage v = |v| exp(i θ).
ExaPF.PowerSystem.VoltageAngle — TypeVoltageAngle <: AbstractNetworkValuesAngle θ of the voltage v = |v| exp(i θ).
ExaPF.PowerSystem.ActivePower — TypeActivePower <: AbstractNetworkValuesActive power P of the complex power S = P + iQ.
ExaPF.PowerSystem.ReactivePower — TypeReactivePower <: AbstractNetworkValuesReactive power Q of the complex power S = P + iQ.
ExaPF.PowerSystem.ActiveLoad — TypeActiveLoad <: AbstractNetworkValuesActive load Pd at buses.
ExaPF.PowerSystem.ReactiveLoad — TypeReactiveLoad <: AbstractNetworkValuesReactive load Qd at buses.
Function to get the range of a given value:
ExaPF.PowerSystem.bounds — Functionbounds(pf::AbstractPowerSystem, attr::AbstractNetworkAttribute, val::AbstractNetworkValues)Return lower and upper bounds corresponding to the admissible values of the AbstractNetworkAttribute attr.
Examples
p_min, p_max = bounds(pf, Generator(), ActivePower())
+v_min, v_max = bounds(pf, Buses(), VoltageMagnitude())
+Settings
This document was generated with Documenter.jl version 1.0.1 on Friday 22 September 2023. Using Julia version 1.9.3.
Given a set of equations F(x) = 0, the Newton-Raphson algorithm for solving nonlinear equations (see below) requires the Jacobian J = jacobian(x) of F. At each iteration a new step dx is computed by solving a linear system. In our case J is sparse and indefinite.
go = true
+ while(go)
+ dx .= jacobian(x)\f(x)
+ x .= x .- dx
+ go = norm(f(x)) < tol ? true : false
+ endThere are two modes of differentiation called forward/tangent or reverse/adjoint. The latter is known in machine learning as backpropagation. The forward mode generates Jacobian-vector product code tgt(x,d) = J(x) * d, while the adjoint mode generates code for the transposed Jacobian-vector product adj(x,y) = (J(x)'*y). We recommend the book Evaluating derivatives: principles and techniques of algorithmic differentiation by Griewank and Walther[1] for a more in-depth introduction to automatic differentiation. The computational complexity of both models favors the adjoint mode if the number of outputs of F is much smaller than the number of inputs size(x) >> size(F), like for example the loss functions in machine learning. However, in our case F is a multivariate vector function from $\mathbb{R}^n$ to $\mathbb{R}^n$, where $n$ is the number of buses.
To avoid a complexity of $\mathcal{O}(n) \cdot cost(F)$ by letting the tangent mode run over all Cartesian basis vectors of $\mathbb{R}^n$, we apply the technique of Jacobian coloring to compress the sparse Jacobian J. Running the tangent mode, it allows to compute columns of the Jacobian concurrently, by combining independent columns in one Jacobian-vector evaluation (see in figure above). For sparsity detection we rely on the greedy algorithm implemented by SparseDiffTools.jl.
Given the sparsity pattern, the forward model is applied through the package ForwardDiff.jl. Given the number of Jacobian colors $c$ we can build our dual type t1s with c directions:
t1s{N} = ForwardDiff.Dual{Nothing,Float64, N} where N}Note that a second-order type t2s can be created naturally by applying the same logic to t1s:
t2s{M,N} = ForwardDiff.Dual{Nothing,t1s{N}, M} where M, N}Finally, this dual type can be ported to both vector types Vector and CuVector:
VT = Vector{Float64}
+VT = Vector{t1s{N}}}
+VT = CuVector{t1s{N}}}Setting VT to either of the three types allows us to instantiate code that has been written using the broadcast operator .
x .= a .* bor accessed in kernels written for KernelAbstractions.jl like for example the power flow equations (here in polar form):
@kernel function residual_kernel!(F, v_m, v_a,
+ ybus_re_nzval, ybus_re_colptr, ybus_re_rowval,
+ ybus_im_nzval, ybus_im_colptr, ybus_im_rowval,
+ pinj, qinj, pv, pq, nbus)
+
+ npv = size(pv, 1)
+ npq = size(pq, 1)
+
+ i = @index(Global, Linear)
+ # REAL PV: 1:npv
+ # REAL PQ: (npv+1:npv+npq)
+ # IMAG PQ: (npv+npq+1:npv+2npq)
+ fr = (i <= npv) ? pv[i] : pq[i - npv]
+ F[i] -= pinj[fr]
+ if i > npv
+ F[i + npq] -= qinj[fr]
+ end
+ @inbounds for c in ybus_re_colptr[fr]:ybus_re_colptr[fr+1]-1
+ to = ybus_re_rowval[c]
+ aij = v_a[fr] - v_a[to]
+ coef_cos = v_m[fr]*v_m[to]*ybus_re_nzval[c]
+ coef_sin = v_m[fr]*v_m[to]*ybus_im_nzval[c]
+ cos_val = cos(aij)
+ sin_val = sin(aij)
+ F[i] += coef_cos * cos_val + coef_sin * sin_val
+ if i > npv
+ F[npq + i] += coef_cos * sin_val - coef_sin * cos_val
+ end
+ end
+endThese two abstractions are a powerful tool that allow us to implement the forward mode in vectorized form where the number of directions or tangent components of a tangent variable are the number of Jacobian colors. We illustrate this in the figure below with a point-wise vector product x .* y
This natural way of computing the compressed Jacobian yields a very high performing code that is portable to any vector architecture, given that a similar package like CUDA.jl exists. We note that similar packages for the Intel Compute Engine and AMD ROCm are in development called oneAPI.jl and AMDGPU.jl, respectively. We expect our package to be portable to AMD and Intel GPUs in the future.
Settings
This document was generated with Documenter.jl version 1.0.1 on Friday 22 September 2023. Using Julia version 1.9.3.
For the purpose of performance regression testing, ExaPF provides a lightweight benchmark script. It allows to test the various configurations for the linear solvers used in the Newton-Raphson algorithm, and run them on a specific hardware. The main julia script benchmark/benchmarks.jl takes all its options from the command line. The benchmark script takes as input a linear solver (e.g. KrylovBICGSTAB), a target architecture as a KernelAbstractions object (CPU or CUDABackend), and a case filename which is included in the ExaData artifact. An exhaustive list of all available linear solvers can be obtained via ExaPF.LinearSolvers.list_solvers.
Running
julia --project benchmark/benchmarks.jl KrylovBICGSTAB CUDABackend case300.myields
KrylovBICGSTAB, CUDABackend, case300.m, 69.0, 3.57, 43.7, trueThe first three fields are the settings of the benchmark run. They are followed by three timings in milliseconds:
To acquire these timings the code is run three times to avoid any precompilation effects. The last field confirms the Newton-Raphson convergence. In case more verbose output is desired, one has to manually set the verbosity in benchmark/benchmarks.jl by changing
powerflow_solver = NewtonRaphson(tol=ntol)to one of the following options:
powerflow_solver = NewtonRaphson(tol=ntol, verbose=VERBOSE_LEVEL_NONE)
+powerflow_solver = NewtonRaphson(tol=ntol, verbose=VERBOSE_LEVEL_LOW)
+powerflow_solver = NewtonRaphson(tol=ntol, verbose=VERBOSE_LEVEL_MEDIUM)
+powerflow_solver = NewtonRaphson(tol=ntol, verbose=VERBOSE_LEVEL_HIGH)A shell script benchmark/benchmarks.sh is provided to gather timings with various canonical configurations and storing them in a file cpu_REV.log and gpu_REF.log, where REV is the sha1 hash of the current checked out ExaPF version.
Settings
This document was generated with Documenter.jl version 1.0.1 on Friday 22 September 2023. Using Julia version 1.9.3.
ExaPF's formalism is based on a vectorized formulation of the power flow problem, as introduced in Lee, Turitsyn, Molzahn, Roald (2020). Throughout this page, we will refer to this formulation as LTMR2020. It is equivalent to the classical polar formulation of the OPF.
In what follows, we denote by $v \in \mathbb{R}^{n_b}$ the voltage magnitudes, $\theta \in \mathbb{R}^{n_b}$ the voltage angles and $p_g, q_g \in \mathbb{R}^{n_g}$ the active and reactive power generations. The active and reactive loads are denoted respectively by $p_d, q_d \in \mathbb{R}^{n_b}$.
The idea is to factorize all nonlinearities inside a basis function depending both on the voltage magnitudes $v$ and voltage angles $\theta$, such that $\psi: \mathbb{R}^{n_b} \times \mathbb{R}^{n_b} \to \mathbb{R}^{2n_\ell + n_b}$. If we introduce the intermediate expressions
\[ \psi_\ell^C(v, \theta) = v^f v^t \cos(\theta_f - \theta_t) \quad \forall \ell = 1, \cdots, n_\ell \\ + \psi_\ell^S(v, \theta) = v^f v^t \sin(\theta_f - \theta_t) \quad \forall \ell = 1, \cdots, n_\ell \\ + \psi_k(v, \theta) = v_k^2 \quad \forall k = 1, \cdots, n_b\]
the basis $\psi$ is defined as
\[ \psi(v, \theta) = [\psi_\ell^C(v, \theta)^\top ~ \psi_\ell^S(v, \theta)^\top ~ \psi_k(v, \theta)^\top ] \, .\]
The basis $\psi$ encodes all the nonlinearities in the problem. For instance, the power flow equations rewrite directly as
\[ \begin{bmatrix} + C_g p_g - p_d \\ + C_g q_g - q_d + \end{bmatrix} + + + \underbrace{ + \begin{bmatrix} + - \hat{G}^c & - \hat{B}^s & -G^d \\ + \hat{B}^c & - \hat{G}^s & B^d + \end{bmatrix} + }_{M} + \psi(v, \theta) + = 0\]
with $C_g \in \mathbb{R}^{n_b \times n_g}$ the bus-generators incidence matrix, and the matrices $B, G$ defined from the admittance matrix $Y_b$ of the network.
Similarly, the line flows rewrite
\[ \begin{bmatrix} + s_p^f \\ s_q^f + \end{bmatrix} + = + \overbrace{ + \begin{bmatrix} + G_{ft} & B_{ft} & G_{ff} C_f^\top \\ + -B_{ft} & G_{ft} & -B_{ff} C_f^\top + \end{bmatrix} + }^{L_{line}^f} + \psi(v, \theta) \\ + \begin{bmatrix} + s_p^t \\ s_q^t + \end{bmatrix} + = + \underbrace{ + \begin{bmatrix} + G_{tf} & B_{tf} & G_{tt} C_t^\top \\ + -B_{tf} & G_{tf} & -B_{tt} C_t^\top + \end{bmatrix} + }_{L_{line}^t} + \psi(v, \theta)\]
with $C_f \in \mathbb{R}^{n_b \times n_\ell}$ the bus-from incidence matrix and $C_t \in \mathbb{R}^{n_b \times n_\ell}$ the bus-to incidence matrix. Then, the line flows constraints write directly with the quadratic expressions:
\[ (s_p^f)^2 + (s_q^f)^2 \leq (s^{max})^2 \quad \, , + (s_p^t)^2 + (s_q^t)^2 \leq (s^{max})^2 \quad \, .\]
Implementing the model LTMR2020 is not difficult once the basis function $\psi$ has been defined. Indeed, if we select a subset of the power flow equations (as usual, associated to the active injections at PV nodes, and active and reactive injections at PQ nodes), we get
\[ C_{eq} p_g + M_{eq} \psi + \tau = 0\]
with $C_{eq}$ defined from the bus-generator incidence matrix $C_g$, $M_{eq}$ a subset of the matrix $M$, $\tau$ a constant depending on the loads in the problem. Note that $C_{eq}$ and $M_{eq}$ are sparse matrices, so the expression can be implemented efficiently with sparse linear algebra operations (2 SpMV operations, 2 vector additions). The same holds true for the line flow constraints, evaluated with 2 SpMV operations:
\[ s^f = L_{line}^f \psi \, , \quad + s^t = L_{line}^t \psi \, .\]
In ExaPF, all nonlinear expressions are written as linear operations depending on the nonlinear basis $\psi$. By doing so, all the unstructured sparsity of the power flow problem is directly handled inside the sparse linear algebra library (cusparse on CUDA GPU, SuiteSparse on the CPU).
In what follows, we detail the implementation of the LTMR2020 model in ExaPF.
We have implemented the LTMR2020 model in ExaPF, both on the CPU and on CUDA GPU. All the operations have been rewritten in a vectorized fashion. Every model depends on inputs we propagate forward with functions. In ExaPF, the inputs will be specified in a NetworkStack <: AbstractStack. The functions will be implemented as AutoDiff.AbstractExpression.
NetworkStackOur three inputs are $(v, \theta, p_g) \in \mathbb{R}^{2n_b + n_g}$ (voltage magnitude, voltage angle, power generations). The basis $\psi$ is considered as an intermediate expression.
We store all inputs in a NetworkStack structure:
struct NetworkStack{VT} <: AbstractStack
+ input::VT
+ vmag::VT # voltage magnitudes (view)
+ vang::VT # voltage angles (view)
+ pgen::VT # active power generations (view)
+ ψ::VT # nonlinear basis ψ(vmag, vang)
+endAll the inputs are specified in the vector input. The three vectors vmag, vang and pgen are views porting on input, and are defined mostly for convenience. By convention the vector input is ordered as [vmag; vang; pgen]:
# Define dimension of the problem
+julia> nbus, ngen, nlines = 3, 2, 4;
+
+julia> stack = ExaPF.NetworkStack(nbus, ngen, nlines, 1, Vector{Float64}, Vector{Float64})
+8-elements NetworkStack{Vector{Float64}}
+
+julia> stack.input
+8-element Vector{Float64}:
+ 0.0
+ 0.0
+ 0.0
+ 0.0
+ 0.0
+ 0.0
+ 0.0
+ 0.0
+
+julia> stack.vmag .= 1;
+
+julia> stack.vang .= 2;
+
+julia> stack.pgen .= 3;
+
+julia> stack.input
+8-element Vector{Float64}:
+ 1.0
+ 1.0
+ 1.0
+ 2.0
+ 2.0
+ 2.0
+ 3.0
+ 3.0The basis vector ψ is an intermediate expression, whose values depend on the inputs.
In the reduced space method, we have to split the variables in a state $x$ and a control $u$. By default, we define
\[ x = (\theta_{pv}, \theta_{pq}, v_{pq}) \, , \quad + x = (v_{ref}, v_{pv}, p_{g,genpv}) \,.\]
and the control, and was not flexible. In the new implementation, we define the state and the control as two mappings porting on the vector stack.input (which itself stores all the inputs in the problem):
julia> nbus, ngen, nlines = 4, 3, 4;
+
+julia> stack = ExaPF.NetworkStack(nbus, ngen, nlines, 1, Vector{Float64}, Vector{Float64});
+
+julia> stack.input .= 1:length(stack.input); # index array input
+
+julia> ref, pv, pq, genpv = [1], [2], [3, 4], [2, 3];
+
+julia> mapx = [nbus .+ pv; nbus .+ pq; pq];
+
+julia> mapu = [ref; pv; 2*nbus .+ genpv];
+
+julia> x = @view stack.input[mapx]
+5-element view(::Vector{Float64}, [6, 7, 8, 3, 4]) with eltype Float64:
+ 6.0
+ 7.0
+ 8.0
+ 3.0
+ 4.0
+
+julia> u = @view stack.input[mapu]
+4-element view(::Vector{Float64}, [1, 2, 10, 11]) with eltype Float64:
+ 1.0
+ 2.0
+ 10.0
+ 11.0By doing so, the values of the state and the control are directly stored inside the NetworkStack structure, avoiding to duplicate values in the memory.
ExaPF implements the different functions required to implement the optimal power flow problem with the polar formulation:
PowerFlowBalance: power flow balance equationsPowerGenerationBounds: bounds on the power generationLineFlows: line flow constraintsCostFunction: quadratic cost functionEach function follows the LTMR2020 model and depends on the basis function $\psi(v, \theta)$, here implemented in the PolarBasis function.
We demonstrate how to use the different functions on the case9 instance. The procedure remains the same for all power network.
julia> polar = ExaPF.load_polar("case9.m");
+
+julia> stack = ExaPF.NetworkStack(polar);
+All the code presented below is agnostic with regards to the specific device (CPU, CUDABackend...) we are using. By default, ExaPF computes the expressions on the CPU. Deporting the computation on a CUDABackend simply translates to instantiate the PolarForm structure on the GPU: polar = PolarForm("case9.m", CUDABackend()).
All functions are following AutoDiff.AbstractExpression's interface. The structure of the network is specified by the PolarForm we pass as an argument in the constructor. For instance, we build a new PolarBasis expression associated to case9 directly as
julia> basis = ExaPF.PolarBasis(polar)
+PolarBasis (AbstractExpression)
+Each expression as a given dimension, given by
julia> length(basis)
+27
+In ExaPF, the inputs and the parameters are stored inside a NetworkStack structure. Evaluating the basis $\psi$ naturally translates to
julia> basis(stack)
+27-element Vector{Float64}:
+ 1.0
+ 1.0
+ 1.0
+ 1.0
+ 1.0
+ 1.0
+ 1.0
+ 1.0
+ 1.0
+ 0.0
+ ⋮
+ 1.0
+ 1.0
+ 1.0
+ 1.0
+ 1.0
+ 1.0
+ 1.0
+ 1.0
+ 1.0This function call allocates a vector psi with 27 elements and evaluates the basis associated to the LTMR2020 model. To avoid unnecessary allocation, one can preallocate the vector psi:
julia> psi = zeros(length(basis)) ;
+
+julia> basis(psi, stack);
+In the LTMR2020 model, the polar basis $\psi(v, \theta)$ depends only on the voltage magnitudes and the voltage angles. However, this is not the case for the other functions (power flow balance, line flows, ...), which all depends on the basis $\psi(v, \theta)$.
In ExaPF, one has to build manually the vectorized expression tree associated to the power flow model. Luckily, evaluating the LTMR2020 simply amounts to compose functions together with the polar basis $\psi(v, \theta)$. ExaPF overloads the function ∘ to compose functions with a PolarBasis instance. The power flow balance can be evaluated as
julia> pflow = ExaPF.PowerFlowBalance(polar) ∘ basis;
+which returns a ComposedExpressions structure.
The function pflow follows the same API, as any regular AutoDiff.AbstractExpression.
julia> n_balance = length(pflow)
+14
+
+julia> round.(pflow(stack); digits=6) # evaluate the power flow balance
+14-element Vector{Float64}:
+ -1.63
+ -0.85
+ 0.0
+ 0.9
+ 0.0
+ 1.0
+ 0.0
+ 1.25
+ -0.167
+ 0.042
+ -0.2835
+ 0.171
+ -0.2275
+ 0.259
+When we evaluate a ComposedExpressions, ExaPF first computes the basis $\psi(v, \theta)$ inside stack.psi, and then ExaPF uses the values in stack.psi to evaluate the final result.
The procedure remains the same if one wants to evaluate the LineFlows or the PowerGenerationBounds. For instance, evaluating the line flows amounts to
julia> line_flows = ExaPF.LineFlows(polar) ∘ basis;
+
+julia> line_flows(stack)
+18-element Vector{Float64}:
+ 0.0
+ 0.006241000000000099
+ 0.0320410000000001
+ 0.0
+ 0.010920249999999961
+ 0.005550250000000068
+ 0.0
+ 0.02340899999999987
+ 0.007743999999999858
+ 0.0
+ 0.006241000000000099
+ 0.0320410000000001
+ 0.0
+ 0.010920249999999961
+ 0.005550250000000068
+ 0.0
+ 0.02340899999999987
+ 0.007743999999999858
+Settings
This document was generated with Documenter.jl version 1.0.1 on Friday 22 September 2023. Using Julia version 1.9.3.
As mentioned before, a linear solver is required to compute the Newton step in
dx .= jacobian(x)\f(x)Our package supports the following linear solvers:
cuSOLVER with csrlsvqr (GPU),Krylov.jl with dqgmres and bicgstab (CPU/GPU),\ operator (CPU),LinearAlgebra.Using only an iterative solver leads to divergence and bad performance due to ill-conditioning of the Jacobian. This is a known phenomenon in power systems. That's why this package comes with a block Jacobi preconditioner that is tailored towards GPUs and is proven to work well with power flow problems.
The Jacobian is partitioned into a dense block diagonal structure using Metis.jl, where each block is inverted to build our preconditioner P.
Compared to incomplete Cholesky and incomplete LU this preconditioner is easily portable to the GPU if the number of blocks is high enough. ExaPF.jl uses the batch BLAS calls from cuBLAS to invert the single blocks.
CUDA.@sync pivot, info = CUDA.CUBLAS.getrf_batched!(blocks, true)
+CUDA.@sync pivot, info, p.cuJs = CUDA.CUBLAS.getri_batched(blocks, pivot)Assuming that other vendors will provide such batched BLAS APIs, this code is portable to other GPU architectures.
Settings
This document was generated with Documenter.jl version 1.0.1 on Friday 22 September 2023. Using Julia version 1.9.3.
The main goal of ExaPF.jl is the solution of optimization problems for electrical power systems in the steady state. The first step in this process is the creation of an object that describes the physics and topology of the power system which ultimately will be mapped into an abstract mathematical optimization problem. In this section we briefly review the power system in the steady state and describe the tools to create and examine power systems in ExaPF.jl.
We usually load the PowerSystem system submodule with the alias PS:
julia> PS = ExaPF.PowerSystem
+The electrical power system is represented as a linear, lumped network which has to satisfy the Kirchhoff laws:
\[ \bm{i} = \bm{Y}\bm{v} \,,\]
where $\bm{i}, \bm{v} \in \mathbb{C}^{N_B}$ are the current and voltage vectors associated to the system and $\bm{Y} \in \mathbb{C}^{N_B \times N_B}$ is the admittance matrix. These equations are often rewritten in terms of apparent powers:
\[ \bm{s} = \bm{p} + j\bm{q} = \textit{diag}(\bm{v^*}) \bm{Y}\bm{v}\]
Using polar representation of the voltage vector, such as $\bm{v} = |v|e^{j \theta}$, each bus $i=1, \cdots, N_B$ must satisfy the power balance equations:
\[\begin{aligned} + p_i &= v_i \sum_{j}^{n} v_j (g_{ij}\cos{(\theta_i - \theta_j)} + b_{ij}\sin{(\theta_i - \theta_j})) \,, \\ + q_i &= v_i \sum_{j}^{n} v_j (g_{ij}\sin{(\theta_i - \theta_j)} - b_{ij}\cos{(\theta_i - \theta_j})) \,. +\end{aligned}\]
where each bus $i$ has variables $p_i, q_i, v_i, \theta_i$ and the topology of the network is defined by a non-negative value of the admittance between two buses $i$ and $j$, $y_{ij} = g_{ij} + ib_{ij}$.
Currently we can create a PS.PowerNetwork object by parsing a MATPOWER data file.
julia> datafile = "case9.m";
+
+julia> ps = PS.load_case(datafile)
+PowerNetwork object with:
+ Buses: 9 (Slack: 1. PV: 2. PQ: 6)
+ Generators: 3.
+Then, using multiple dispatch, we have defined a set of abstract data types and getter functions which allow us to retrieve information from the PowerNetwork object
julia> PS.get(ps, PS.NumberOfPQBuses())
+6
+
+julia> PS.get(ps, PS.NumberOfPVBuses())
+2
+
+julia> PS.get(ps, PS.NumberOfSlackBuses())
+1Settings
This document was generated with Documenter.jl version 1.0.1 on Friday 22 September 2023. Using Julia version 1.9.3.
This page introduces the first steps to set up ExaPF.jl. We show how to load a power network instance and how to solve the power flow equations both on the CPU and on the GPU. The full script is implemented in test/quickstart.jl.
We start by importing CUDA and KernelAbstractions:
using CUDA
+using KernelAbstractionsThen, we load ExaPF and its submodules with
using ExaPF
+import ExaPF: AutoDiff
+const PS = ExaPF.PowerSystem
+const LS = ExaPF.LinearSolversExaPF loads instances from the pglib-opf benchmark. ExaPF contains an artifact defined in Artifacts.toml that is built from the ExaData repository containing Exascale Computing Project relevant test cases. You may set a data file using
datafile = joinpath(artifact"ExaData", "ExaData", "case1354.m")The powerflow equations can be solved in three lines of code, as
julia> polar = ExaPF.PolarForm(datafile, CPU()) # Load dataPolar formulation (instantiated on device CPU(false)) +Network characteristics: + #buses: 1354 (#slack: 1 #PV: 259 #PQ: 1094) + #generators: 260 + #lines: 1991 +giving a mathematical formulation with: + #controls: 519 + #states : 2447julia> stack = ExaPF.NetworkStack(polar) # Load variables2968-elements NetworkStack{Vector{Float64}}julia> convergence = run_pf(polar, stack; verbose=1)#it 0: 1.15103e+02 +#it 1: 1.50328e+01 +#it 2: 5.88242e-01 +#it 3: 4.88493e-03 +#it 4: 1.39924e-06 +#it 5: 1.52178e-11 +Power flow has converged: true + * #iterations: 5 + * Time Jacobian (s) ........: 0.4035 + * Time linear solver (s) ...: 0.0207 + * Time total (s) ...........: 0.8188
Implicitly, ExaPF has just proceed to the following operations:
stackThis compact syntax allows to solve quickly any powerflow equations in a few lines a code. However, in most case, the user may want more coarse grained control on the different objects manipulated.
In what follows, we detail step by step the detailed procedure to solve the powerflow equations.
We start by importing a MATPOWER instance to a ExaPF.PowerSystem.PowerNetwork object:
julia> pf = PS.PowerNetwork(datafile)PowerNetwork object with: + Buses: 1354 (Slack: 1. PV: 259. PQ: 1094) + Generators: 260.
The different fields of the object pf specify the characteristics of the network. For instance, we can retrieve the number of buses or get the indexing of the PV buses with
julia> nbus = PS.get(pf, PS.NumberOfBuses())1354julia> pv_indexes = pf.pv;
However, a ExaPF.PowerSystem.PowerNetwork object stores only the physical attributes of the network. To choose a given mathematical formulation, we need to pass the object pf to an ExaPF.AbstractFormulation layer. Currently, only the polar formulation is provided with the ExaPF.PolarForm structure. In the future, other formulations (e.g. RectangularForm) may be implemented as well.
To solve the powerflow equations, we need to choose a given mathematical formulation for the equations of the network. To each formulation corresponds a given state $x$ and control $u$. Using polar representation of the voltage vector, such as $\bm{v} = |v|e^{j \theta}$, each bus $i=1, \cdots, N_B$ must satisfy the power balance equations:
\[\begin{aligned} + p_i &= v_i \sum_{j}^{n} v_j (g_{ij}\cos{(\theta_i - \theta_j)} + b_{ij}\sin{(\theta_i - \theta_j})) \,, \\ + q_i &= v_i \sum_{j}^{n} v_j (g_{ij}\sin{(\theta_i - \theta_j)} - b_{ij}\cos{(\theta_i - \theta_j})) \,. +\end{aligned}\]
The powerflow equations rewrite in the abstract mathematical formalism:
\[g(x, u) = 0.\]
For a given control $u$, solving the powerflow equations resumes to find a state $x(u)$ such that $g(x(u), u) = 0$.
To this goal, ExaPF.jl implements a Newton-Raphson algorithm that allows to solve the powerflow equations in a few lines of code. We first instantiate a PolarForm object to adopt a polar formulation as a model:
julia> polar = ExaPF.PolarForm(pf, CPU())Polar formulation (instantiated on device CPU(false)) +Network characteristics: + #buses: 1354 (#slack: 1 #PV: 259 #PQ: 1094) + #generators: 260 + #lines: 1991 +giving a mathematical formulation with: + #controls: 519 + #states : 2447
Note that the constructor ExaPF.PolarForm takes as input a ExaPF.PowerSystem.PowerNetwork object and a KernelAbstractions.jl device (here set to CPU() by default). We will explain in the next section how to load a ExaPF.PolarForm object on the GPU with the help of a CUDABackend().
The Newton-Raphson solves the equation $g(x, u) = 0$ in an iterative fashion. The algorithm solves at each step the linear equation:
\[ x_{k+1} = - (\nabla_x g_k)^{-1} g(x_k, u).\]
Hence, the algorithm requires the following elements:
The variable $x$ is instantiated as:
julia> stack = ExaPF.NetworkStack(polar)2968-elements NetworkStack{Vector{Float64}}
The function $g$ is implemented using ExaPF's custom modeler:
julia> basis = ExaPF.PolarBasis(polar)PolarBasis (AbstractExpression)julia> powerflow = ExaPF.PowerFlowBalance(polar) ∘ basisExaPF.ComposedExpressions{ExaPF.PolarBasis{Vector{Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}}, ExaPF.PowerFlowBalance{Vector{Float64}, SparseArrays.SparseMatrixCSC{Float64, Int64}}}(PolarBasis (AbstractExpression), PowerFlowBalance (AbstractExpression))
The Jacobian $\nabla_x g$ is evaluated automatically using forward-mode AutoDiff:
julia> mapx = ExaPF.mapping(polar, State());julia> jx = ExaPF.Jacobian(polar, powerflow, mapx)A AutoDiff Jacobian for ExaPF.ComposedExpressions{ExaPF.PolarBasis{Vector{Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}}, ExaPF.PowerFlowBalance{Vector{Float64}, SparseArrays.SparseMatrixCSC{Float64, Int64}}}(PolarBasis (AbstractExpression), PowerFlowBalance (AbstractExpression)) +Number of Jacobian colors: 25
The (direct) linear solver can be instantiated directly as
julia> linear_solver = LS.DirectSolver(jx.J);
Let's explain further these three objects.
stack is a AbstractStack storing all the variables attached to the formulation polar::PolarForm.jx is a Jacobian structure which allows the solver to compute efficiently the Jacobian of the powerflow equations $\nabla_x g$ using AutoDiff.linear_solver specifies the linear algorithm uses to solve the linear system $(\nabla_x g_k) x_{k+1} = g(x_k, u)$. By default, we use direct linear algebra.In the AutoDiff Jacobian jx, the evaluation of the Jacobian $J$ is stored in jx.J:
julia> jac = jx.J2447×2447 SparseArrays.SparseMatrixCSC{Float64, Int64} with 15803 stored entries: +⎡⢿⣷⣲⣷⣯⣷⣣⢯⣮⣌⢫⡿⣿⣲⢡⣾⡾⣶⣰⣷⡧⣾⣵⣿⢺⣵⣧⣙⣿⢿⣟⡮⢴⣗⣷⣶⣾⣿⡱⡧⎤ +⎢⢼⣾⣿⣿⣿⡺⣶⣿⣟⢷⣲⣿⣷⣿⣳⡟⢹⡿⣿⢿⢓⣻⡗⣾⢼⡿⡿⣆⣾⣾⣯⣗⢿⠃⡻⣿⡿⣟⢚⣧⎥ +⎢⢯⣿⣻⡻⢟⣵⢿⣟⣇⣾⣵⣿⡿⡷⣿⣿⢿⣿⢿⢟⣿⡽⢿⣾⢟⣾⣷⣫⣾⣿⣿⡾⣯⣾⣟⡿⡷⣿⣿⠿⎥ +⎢⡭⣞⣼⣿⣿⢷⡿⣯⣧⣿⢿⣾⣷⣷⡫⣧⣻⣻⣵⢾⡹⣻⡿⣿⢿⣿⣹⣾⣿⣿⣶⡟⣻⣌⣿⣮⣷⡏⣿⡇⎥ +⎢⡊⢿⢿⣝⣩⣽⣭⣿⣿⢟⣻⣿⡷⢿⢽⣯⣿⢾⣼⡻⢽⣼⢯⡻⣽⢿⡿⣟⣿⣿⠿⣿⣿⣽⡷⣿⣽⡯⣧⣿⎥ +⎢⣯⡶⣼⣾⣵⣿⣻⣷⣿⣾⡛⣬⣛⣿⢻⣽⠽⣯⣿⣳⣙⣿⣾⣝⣿⣴⡿⣿⢳⣽⣻⡯⣯⡿⢧⣻⣾⣏⣻⢗⎥ +⎢⢻⣻⣽⣿⢿⡯⢽⣿⣽⣏⣿⣼⡿⣯⣞⣷⡿⣹⡽⣗⣾⣟⢵⡭⣿⡯⣷⣹⣧⣻⢿⣷⣳⣳⠯⣿⣯⣷⣾⡝⎥ +⎢⣡⣶⣽⠾⣿⣿⠯⣮⡷⣷⣟⣶⢾⣽⣻⣾⣷⣷⢿⣿⣾⣳⣿⣿⢵⣾⢿⣾⣻⣦⢯⣞⣿⣿⣿⢿⣟⣿⣟⣜⎥ +⎢⢺⣯⣷⡶⣿⣷⣿⣺⣻⣟⡷⣧⣟⣫⢽⣿⡻⣮⣿⣳⢻⣿⣷⢿⢟⣞⣯⣷⣻⣾⣓⣩⣿⣟⢿⣾⣟⡮⣺⣿⎥ +⎢⢴⣾⣿⣟⣿⢗⣱⣟⣶⡻⢿⣻⢷⢯⣿⣷⢿⣻⢿⢗⣚⣞⣾⣒⣿⣷⣖⡷⣿⡿⡼⣿⣿⣾⣿⣯⣿⣖⢺⡹⎥ +⎢⣩⣯⣽⣰⣟⡿⣷⣪⣓⣷⣷⣼⣾⢿⢾⣻⣿⣶⣺⢼⡿⣯⣿⢿⣗⣝⣹⣿⢷⣷⡾⡷⢟⢿⣷⣲⡯⣞⣿⣿⎥ +⎢⣵⣿⣹⣭⣻⣷⣿⣯⣯⡳⣞⢿⡕⡷⣿⣿⣽⣟⢺⢻⣿⣟⣿⣿⢯⣟⣟⢶⣾⣻⣻⣿⣷⢯⣿⣟⡟⣿⡽⡻⎥ +⎢⢞⣶⣶⡷⣻⣵⣿⣷⣷⣟⢛⣿⡿⡿⣱⣷⣻⢵⢿⣿⣝⢽⣯⢷⣿⣿⣽⣟⣻⣿⢿⠏⣹⣦⡿⢶⣿⣏⡿⢇⎥ +⎢⣍⢻⠻⢯⡽⣻⣳⣾⣿⢯⣿⣯⣝⣻⣻⣷⢯⣿⢼⡽⣷⣾⢻⣝⣷⢿⣟⣽⡻⣭⣛⣽⢿⡿⡿⣳⣾⣗⣷⡻⎥ +⎢⣿⣟⣺⣿⣾⣿⣿⣿⣿⣿⣝⣶⣭⣻⠻⣾⣻⣾⣿⡿⢽⣷⣾⣻⣿⣾⡟⣮⣻⣾⣝⡷⢷⡟⣧⣾⣿⡯⢿⣲⎥ +⎢⡻⡽⢯⢿⣻⡿⣼⠿⣿⣧⡿⡾⢿⣷⣫⢷⡝⣸⣶⣯⢾⡯⣿⣾⡿⠗⣟⣼⢷⡽⣿⣿⣹⣿⣗⡿⣦⣷⣿⣎⎥ +⎢⢴⢷⠿⠓⣫⣿⡛⢾⣟⣿⣯⡿⢽⣺⣿⣿⣿⢿⣻⣿⣿⣕⡽⣟⠳⣾⣿⡷⣽⠷⣷⣾⣿⣿⣿⢿⣿⣿⣯⣶⎥ +⎢⢹⣿⣿⣮⣿⡽⡻⣿⣽⣯⣭⣳⣯⣧⣿⣟⣻⣷⡿⣿⢹⣻⣿⢿⢻⣏⢿⣫⣩⣿⣽⡽⣿⣟⣵⣿⣿⡷⣽⣷⎥ +⎢⣾⣿⣿⢯⣽⣯⡽⠿⡷⡿⡾⢿⢯⣿⣿⣽⡻⡽⢻⢿⣫⢯⣿⣭⡿⢿⢾⢿⡿⡿⢬⣿⣿⣿⢿⡿⣵⣿⣹⣝⎥ +⎣⠵⡮⠾⣴⣿⡟⠿⠿⣭⣿⢿⢞⣞⠿⣛⢽⣾⣾⣞⡲⣿⣿⣷⡫⠿⢏⣽⡻⢻⣳⡻⢿⢫⣿⢷⣿⣗⢾⣿⣿⎦
This matrix is at the basis of the powerflow algorithm. At each iteration, the AutoDiff backend updates the nonzero values in the sparse Jacobian jx and solve the associated linear system to compute the next descent direction.
The procedure is implemented in the nlsolve! function, which uses a Newton-Raphson algorithm to solve the powerflow equations. The Newton-Raphson algorithm is specified as:
julia> pf_algo = NewtonRaphson(; verbose=1, tol=1e-10)NewtonRaphson(20, 1.0e-10, 1)
Then, we can solve the powerflow equations simply with
julia> convergence = ExaPF.nlsolve!(pf_algo, jx, stack; linear_solver=linear_solver)#it 0: 1.15103e+02 +#it 1: 1.50328e+01 +#it 2: 5.88242e-01 +#it 3: 4.88493e-03 +#it 4: 1.39924e-06 +#it 5: 1.52178e-11 +Power flow has converged: true + * #iterations: 5 + * Time Jacobian (s) ........: 0.0045 + * Time linear solver (s) ...: 0.0221 + * Time total (s) ...........: 0.0269
Here, the algorithm solves the powerflow equations in 5 iterations. The algorithm modifies the values of stack inplace, to avoid any unnecessary memory allocations.
Now, how can we deport the resolution on the GPU? The procedure looks exactly the same. It suffices to initiate a new ExaPF.PolarForm object, but on the GPU:
julia> polar_gpu = ExaPF.PolarForm(pf, CUDABackend())Polar formulation (instantiated on device CUDA.CUDAKernels.CUDABackend(false, false)) +Network characteristics: + #buses: 1354 (#slack: 1 #PV: 259 #PQ: 1094) + #generators: 260 + #lines: 1991 +giving a mathematical formulation with: + #controls: 519 + #states : 2447
polar_gpu will load all the structures it needs on the GPU, to avoid unnecessary movements between the host and the device. We can load the other structures directly on the GPU with:
julia> stack_gpu = ExaPF.NetworkStack(polar_gpu)2968-elements NetworkStack{CUDA.CuArray{Float64, 1, CUDA.Mem.DeviceBuffer}}julia> basis_gpu = ExaPF.PolarBasis(polar_gpu)PolarBasis (AbstractExpression)julia> pflow_gpu = ExaPF.PowerFlowBalance(polar_gpu) ∘ basis_gpuExaPF.ComposedExpressions{ExaPF.PolarBasis{CUDA.CuArray{Int64, 1}, CUDA.CUSPARSE.CuSparseMatrixCSR{Float64, Int32}}, ExaPF.PowerFlowBalance{CUDA.CuArray{Float64, 1}, CUDA.CUSPARSE.CuSparseMatrixCSR{Float64, Int32}}}(PolarBasis (AbstractExpression), PowerFlowBalance (AbstractExpression))julia> jx_gpu = ExaPF.Jacobian(polar_gpu, pflow_gpu, mapx)A AutoDiff Jacobian for ExaPF.ComposedExpressions{ExaPF.PolarBasis{CUDA.CuArray{Int64, 1}, CUDA.CUSPARSE.CuSparseMatrixCSR{Float64, Int32}}, ExaPF.PowerFlowBalance{CUDA.CuArray{Float64, 1}, CUDA.CUSPARSE.CuSparseMatrixCSR{Float64, Int32}}}(PolarBasis (AbstractExpression), PowerFlowBalance (AbstractExpression)) +Number of Jacobian colors: 25julia> linear_solver = LS.DirectSolver(jx_gpu.J)ExaPF.LinearSolvers.DirectSolver{Nothing}(nothing)
Then, solving the powerflow equations on the GPU directly translates as
julia> convergence = ExaPF.nlsolve!(pf_algo, jx_gpu, stack_gpu; linear_solver=linear_solver)#it 0: 1.15103e+02 +#it 1: 1.50328e+01 +#it 2: 5.88242e-01 +#it 3: 4.88493e-03 +#it 4: 1.39924e-06 +#it 5: 9.96622e-12 +Power flow has converged: true + * #iterations: 5 + * Time Jacobian (s) ........: 2.6411 + * Time linear solver (s) ...: 0.3088 + * Time total (s) ...........: 3.7838
Note that we get exactly the same iterations as when we solve the power flow equations on the CPU.
By default, the algorithm runs with a direct solver, which might be inefficient for large problems. To overcome this issue, ExaPF implements a wrapper for different iterative algorithms (GMRES, BICGSTAB).
The performance of iterative solvers is usually improved if we use a preconditioner. ExaPF.jl implements an overlapping Schwarz preconditioner, tailored for GPU usage. To build an instance with 8 blocks, just write
julia> npartitions = 8;julia> jac_gpu = jx_gpu.J;julia> precond = LS.BlockJacobiPreconditioner(jac_gpu, npartitions, CUDABackend());
You can attach the preconditioner to an BICGSTAB algorithm simply as
julia> linear_solver = ExaPF.KrylovBICGSTAB(jac_gpu; P=precond);
(this will use the BICGSTAB algorithm implemented in Krylov.jl).
We need to update accordingly the tolerance of the Newton-Raphson algorithm (the iterative solver is less accurate than the direct solver):
julia> pf_algo = NewtonRaphson(; verbose=1, tol=1e-7)NewtonRaphson(20, 1.0e-7, 1)
We reset the variables to their initial values:
julia> ExaPF.init!(polar_gpu, stack_gpu)
Then, solving the power flow with the iterative solvers directly translates to one call to nlsolve!:
julia> convergence = ExaPF.nlsolve!(pf_algo, jx_gpu, stack_gpu; linear_solver=linear_solver)#it 0: 1.15103e+02 +#it 1: 1.50328e+01 +#it 2: 5.88242e-01 +#it 3: 4.88493e-03 +#it 4: 1.39924e-06 +#it 5: 6.63821e-11 +Power flow has converged: true + * #iterations: 5 + * Time Jacobian (s) ........: 0.0031 + * Time linear solver (s) ...: 6.2590 + * Time total (s) ...........: 6.2635
Settings
This document was generated with Documenter.jl version 1.0.1 on Friday 22 September 2023. Using Julia version 1.9.3.
ExaPF provides a way to evaluate the expressions by blocks, opening the way to introduce more parallelism in the code.
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.20808128302365436
+ 0.380939002512915
+ 0.2086878043432392
+ 0.1460768320108018
+ 0.7285197176703684
+ 0.7700845730585193
+ 0.42044624137648023
+ 0.25391611879510856
+ 0.8807707460953512
+ 0.9843937115059406
+ ⋮
+ 0.1436296331486917
+ 0.2300357962635473
+ 0.14804970047711885
+ 0.6670337145276783
+ 0.04286817446758806
+ 0.23098603656336025
+ 0.4059550496984342
+ 0.5448914927590742
+ 0.2572340418875648The 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.871771 0.685618 0.21964 0.376286 … 0.767582 0.557717 0.652547
+ 0.172378 0.207749 0.265254 0.394895 0.606452 0.233679 0.479319
+ 0.740584 0.426237 0.901151 0.650416 0.526265 0.447807 0.981291
+ 0.868898 0.513159 0.130942 0.650565 0.81149 0.199634 0.583923
+ 0.0884298 0.500398 0.635016 0.180544 0.0896458 0.332904 0.7909
+ 0.7995 0.299343 0.315024 0.922389 … 0.844315 0.00204092 0.081186
+ 0.256418 0.856264 0.252553 0.992181 0.273027 0.930412 0.0158826
+ 0.392837 0.885696 0.529127 0.530639 0.800562 0.336638 0.770553
+ 0.548838 0.477367 0.378557 0.875426 0.605468 0.427606 0.0493593ExaPF 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.45762 -1.42225 -1.36475 … -1.02355 -1.39632 -1.15068
+ -0.109416 -0.423763 0.0511506 -0.323735 -0.402193 0.131291
+ 0.868898 0.513159 0.130942 0.81149 0.199634 0.583923
+ 0.0884298 0.500398 0.635016 0.0896458 0.332904 0.7909
+ 0.7995 0.299343 0.315024 0.844315 0.00204092 0.081186
+ 0.256418 0.856264 0.252553 … 0.273027 0.930412 0.0158826
+ 0.392837 0.885696 0.529127 0.800562 0.336638 0.770553
+ 0.548838 0.477367 0.378557 0.605468 0.427606 0.0493593
+ -0.0209232 0.318834 0.173779 -0.124165 0.443083 0.500034
+ 0.47052 0.600145 -0.250967 -0.125072 0.624712 -0.215132
+ 0.486585 -0.141536 0.368738 … 0.320649 0.423069 -0.052514
+ 0.241446 0.0204472 0.24379 0.427511 -0.151737 0.226955
+ 0.0264161 0.267377 -0.011641 0.0250835 0.033149 0.317391
+ 0.639771 0.519448 0.230413 0.151449 -0.22263 0.016234We get $N$ different results for the power flow balance equations, depending on which scenario we are on.
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: 5
+ * Time Jacobian (s) ........: 0.0002
+ * Time linear solver (s) ...: 0.0273
+ * Time total (s) ...........: 0.0277
+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.949149 0.937853 0.98378 0.957957 0.992592 0.958013 0.975337
+ 0.913254 0.899809 0.986506 0.934134 0.989741 0.912533 0.981569
+ 0.945667 0.972952 0.973366 0.957304 … 0.968594 0.966219 0.98814
+ 0.939979 0.957841 0.964845 0.911051 0.950752 0.979927 0.969126
+ 0.960926 0.958658 0.981665 0.945627 0.974497 0.98877 0.974117
+ 0.905152 0.906513 0.964963 0.936795 0.972029 0.975694 0.972684When 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: 5
+ * Time Jacobian (s) ........: 0.0014
+ * Time linear solver (s) ...: 0.0569
+ * Time total (s) ...........: 0.0594
+Settings
This document was generated with Documenter.jl version 1.0.1 on Friday 22 September 2023. Using Julia version 1.9.3.
ExaPF implements a power flow solver in the function run_pf. Under the hood, the function run_pf calls the function nlsolve! which uses a Newton-Raphson algorithm to solve iteratively the system of nonlinear equations
\[g(x, p) = 0\]
where $g: \mathbb{R}^{n_x} \times \mathbb{R}^{n_p} \to \mathbb{R}^{n_x}$ is a nonlinear function encoding the power flow equations.
At a fixed $p$, solving the power flow amounts to find a state $x$ such that $g(x, p) = 0$. At iteration $$k$, the Newton-Raphson algorithm finds the next iterate by solving the linear system
\[(\nabla_x g_k) \Delta x = - g_k\]
and setting $x_{k+1} = x_{k} + \Delta x_k$. The Jacobian $\nabla_x g_k = \nabla_x (x_k, p)$ is computed automatically in sparse format using AutoDiff.
Hence, solving the power flow equations amounts to solve a sequence of sparse linear systems. When a direct solver is employed, the system is solved in two steps. First, a LU factorization of the matrix $\nabla_x g$ is computed: we find a lower and an upper triangular matrices $L$ and $U$ as well as two permutation matrices $P$ and $Q$ such that
\[P (\nabla_x g) Q = LU\]
Once the matrix factorized, solving the linear system just translates to perform two backsolves with the triangular matrices $L$ and $U$.
This method is usually efficient, as the power flow Jacobian is super sparse (less than 1% of nonzeroes) and its sparsity pattern is fixed, so we have to compute the symbolic factorization of the system only once.
By default, ExaPF employs the linear solver UMFPACK to solve the linear system, as UMFPACK is shipped automatically in Julia.
In the LinearSolvers submodule, this is how the wrapper is implemented:
struct DirectSolver{Fac} <: AbstractLinearSolver
+ factorization::Fac
+end
+DirectSolver(J::AbstractMatrix) = DirectSolver(lu(J))
+By default, the constructor takes as input the initial Jacobian J and factorizes it by calling lu(J), which in Julia translates to a factorization with UMFPACK. Then, inside the function nlsolve! we refactorize the matrix at each iteration by calling the function LinearSolvers.update!
function update!(s::DirectSolver, J::AbstractMatrix)
+ LinearAlgebra.lu!(s.factorization, J) # Update factorization inplace
+endThis function uses the function LinearAlgebra.lu! to update the factorization inplace. The backsolve is computed by calling the LinearAlgebra.ldiv! function:
function ldiv!(s::DirectSolver, y::AbstractVector, J::AbstractMatrix, x::AbstractVector)
+ LinearAlgebra.ldiv!(y, s.factorization, x) # Forward-backward solve
+ return 0
+endWe notice that the code has been designed to support any factorization routines implementing the two routines LinearAlgebra.lu! and LinearAlgebra.ldiv!.
Before comparing with other linear solvers, we solve a large scale power flow instance with UMFPACK to give us a reference.
polar = ExaPF.load_polar("case9241pegase.m")
+stack = ExaPF.NetworkStack(polar)
+pf_solver = NewtonRaphson(tol=1e-10, verbose=2) # power flow solver
+func = ExaPF.PowerFlowBalance(polar) ∘ ExaPF.PolarBasis(polar) # power flow func
+jx = ExaPF.Jacobian(polar, func, State()) # init AD
+ExaPF.nlsolve!(pf_solver, jx, stack)Power flow has converged: true
+ * #iterations: 7
+ * Time Jacobian (s) ........: 1.1377
+ * Time linear solver (s) ...: 0.3192
+ * Time total (s) ...........: 1.4619
+KLU is an efficient sparse linear solver, initially designed for circuit simulation problems. It is often considered as one of the state-of-the-art linear solver to solve power flow problems. Conveniently, KLU is wrapped in Julia with the package KLU.jl. KLU.jl implements a proper interface to use KLU. We just have to implement a forgiving function for LinearAlgebra.lu!
LinearAlgebra.lu!(K::KLU.KLUFactorization, J) = KLU.klu!(K, J)Then, we are ready to solve a power flow with KLU using our current abstraction. One has just to create a new instance of LS.DirectSolver:
klu_factorization = KLU.klu(jx.J)
+klu_solver = LS.DirectSolver(klu_factorization)ExaPF.LinearSolvers.DirectSolver{KLU.KLUFactorization{Float64, Int64}}(KLU.KLUFactorization{Float64, Int64}(KLU.LibKLU.klu_l_common_struct(0.001, 1.2, 1.2, 10.0, 0.0, 1, 0, 2, Ptr{Nothing} @0x0000000000000000, Ptr{Nothing} @0x0000000000000000, 1, 0, 0, 17036, 17036, 17036, 0, -1.0, -1.0, -1.0, -1.0, 0.0, 0x00000000005a0030, 0x0000000000760150), Ptr{Nothing} @0x000000004ce0efe0, Ptr{Nothing} @0x000000004bebf320, 17036, [0, 4, 12, 19, 24, 27, 37, 40, 44, 54 … 129350, 129359, 129363, 129369, 129380, 129384, 129390, 129402, 129408, 129412], [0, 383, 5548, 13344, 1, 458, 4027, 5306, 6117, 11823 … 6009, 6175, 9238, 13805, 13971, 17034, 8879, 9239, 16675, 17035], [1050.4551624894714, -629.5373899630948, -421.251258755022, 88.40062041224645, 995.2128749049234, -479.4136240982223, -138.29482949968113, -145.959631463846, -230.919755492703, 22.972728768573532 … -480.4784394426574, -39.711370520688675, 520.164296001656, -3360.900457400469, -208.11090119999253, 3568.9822586698524, -39.40403545216682, 39.40522422897138, -242.80215386396614, 242.80196093598838]))and pass it to nlsolve!:
ExaPF.init!(polar, stack) # reinit stack
+ExaPF.nlsolve!(pf_solver, jx, stack; linear_solver=klu_solver)Power flow has converged: true
+ * #iterations: 7
+ * Time Jacobian (s) ........: 0.2740
+ * Time linear solver (s) ...: 0.0330
+ * Time total (s) ...........: 0.3120
+We observe KLU reduces considerably the time spent in the linear solver.
cusolverRF is an efficient LU refactorization routine implemented in CUDA. It is wrapped in Julia inside the package CUSOLVERRF.jl:
using CUSOLVERRFThe principle is the following: the initial symbolic factorization is computed on the CPU with the routine chosen by the user. Then, each time we have to refactorize a matrix with the same sparsity pattern, we can recompute the numerical factorization entirely on the GPU. In practice, this solver is efficient at refactorizing a given matrix if the sparsity is significant.
This is of direct relevance for us, as (i) the sparsity of the power flow Jacobian doesn't change along the Newton iterations and (ii) the Jacobian is super-sparse. In ExaPF, it is the linear solver of choice when it comes to solve the power flow entirely on the GPU.
CUSOLVERRF.jl follows the LinearAlgebra's interface, so we can use it directly in ExaPF. We first have to instantiate everything on the GPU:
using CUDA
+polar_gpu = ExaPF.load_polar("case9241pegase.m", CUDABackend())
+stack_gpu = ExaPF.NetworkStack(polar_gpu)
+func_gpu = ExaPF.PowerFlowBalance(polar_gpu) ∘ ExaPF.PolarBasis(polar_gpu)
+jx_gpu = ExaPF.Jacobian(polar_gpu, func_gpu, State()) # init ADA AutoDiff Jacobian for ExaPF.ComposedExpressions{ExaPF.PolarBasis{CUDA.CuArray{Int64, 1}, CUDA.CUSPARSE.CuSparseMatrixCSR{Float64, Int32}}, ExaPF.PowerFlowBalance{CUDA.CuArray{Float64, 1}, CUDA.CUSPARSE.CuSparseMatrixCSR{Float64, Int32}}}(PolarBasis (AbstractExpression), PowerFlowBalance (AbstractExpression))
+Number of Jacobian colors: 85We can instantiate a new cusolverRF's instance as
rf_fac = CUSOLVERRF.RFLU(jx_gpu.J)
+rf_solver = LS.DirectSolver(rf_fac)ExaPF.LinearSolvers.DirectSolver{CUSOLVERRF.RFLU{Float64, CUSOLVERRF.CuSparseSVDeprecated, CUSOLVERRF.CuSparseSMDeprecated}}(CUSOLVERRF.RFLU{Float64, CUSOLVERRF.CuSparseSVDeprecated, CUSOLVERRF.CuSparseSMDeprecated}
+)Then, we are able to solve the power flow entirely on the GPU, simply as
ExaPF.nlsolve!(pf_solver, jx_gpu, stack_gpu; linear_solver=rf_solver)Power flow has converged: true
+ * #iterations: 7
+ * Time Jacobian (s) ........: 1.3818
+ * Time linear solver (s) ...: 0.0822
+ * Time total (s) ...........: 1.4656
+Settings
This document was generated with Documenter.jl version 1.0.1 on Friday 22 September 2023. Using Julia version 1.9.3.