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 0.27.25 on Saturday 2 September 2023. Using Julia version 1.8.5.
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 Thursday 21 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 0.27.25 on Saturday 2 September 2023. Using Julia version 1.8.5.
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 Thursday 21 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
+AutoDiff · ExaPF.jl AutoDiff
Variables
ExaPF.AutoDiff.AbstractStack — TypeAbstractStack{VT}
Abstract variable storing the inputs and the intermediate values in the expression tree.
sourceExpressions
ExaPF.AutoDiff.AbstractExpression — TypeAbstractExpression
Abstract 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.
sourceExaPF.AutoDiff.adjoint! — Functionadjoint!(expr::AbstractExpression, ∂stack::AbstractStack{VT}, stack::AbstractStack{VT}, ̄v::VT) where VT<:AbstractArray
Compute 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!.
sourceFirst and second-order derivatives
ExaPF.AutoDiff.AbstractJacobian — TypeAbstractJacobian
Automatic differentiation for the compressed Jacobian of any nonlinear constraint $h(x)$.
sourceExaPF.AutoDiff.AbstractHessianProd — TypeAbstractHessianProd
Returns the adjoint-Hessian-vector product $λ^⊤ H v$ of any nonlinear constraint $h(x)$.
sourceExaPF.AutoDiff.AbstractFullHessian — TypeAbstractHessianProd
Full sparse Hessian $H$ of any nonlinear constraint $h(x)$.
sourceUtils
ExaPF.AutoDiff.seed! — Functionseed!(
+
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.
sourceExaPF.AutoDiff.adjoint! — Functionadjoint!(expr::AbstractExpression, ∂stack::AbstractStack{VT}, stack::AbstractStack{VT}, ̄v::VT) where VT<:AbstractArray
Compute 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!.
sourceFirst and second-order derivatives
ExaPF.AutoDiff.AbstractJacobian — TypeAbstractJacobian
Automatic differentiation for the compressed Jacobian of any nonlinear constraint $h(x)$.
sourceExaPF.AutoDiff.AbstractHessianProd — TypeAbstractHessianProd
Returns the adjoint-Hessian-vector product $λ^⊤ H v$ of any nonlinear constraint $h(x)$.
sourceExaPF.AutoDiff.AbstractFullHessian — TypeAbstractHessianProd
Full sparse Hessian $H$ of any nonlinear constraint $h(x)$.
sourceUtils
ExaPF.AutoDiff.seed! — Functionseed!(
H::AbstractHessianProd,
v::AbstractVector{T},
-) where {T}
Seed the duals with v to compute the Hessian vector product $λ^⊤ H v$.
sourceExaPF.AutoDiff.seed_coloring! — Functionseed_coloring!(
+) where {T}
Seed the duals with v to compute the Hessian vector product $λ^⊤ H v$.
sourceExaPF.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$.
sourceExaPF.AutoDiff.partials! — Functionpartials!(jac::AbstractJacobian)
Extract partials from Jacobian jac in jac.J.
sourcepartials!(hess::AbstractFullHessian)
Extract partials from Hessian hess into hess.H.
sourceExaPF.AutoDiff.set_value! — Functionset_value!(
+)
Seed the duals with the coloring based seeds to compute the Jacobian or Hessian $M$.
sourceExaPF.AutoDiff.partials! — Functionpartials!(jac::AbstractJacobian)
Extract partials from Jacobian jac in jac.J.
sourcepartials!(hess::AbstractFullHessian)
Extract partials from Hessian hess into hess.H.
sourceExaPF.AutoDiff.set_value! — Functionset_value!(
jac,
primals::AbstractVector{T}
-) where {T}
Set values of ForwardDiff.Dual numbers in jac to primals.
sourceSettings
This document was generated with Documenter.jl version 0.27.25 on Saturday 2 September 2023. Using Julia version 1.8.5.
Set values of ForwardDiff.Dual numbers in jac to primals.
Settings
This document was generated with Documenter.jl version 1.0.1 on Thursday 21 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;
+Polar formulation · ExaPF.jl Polar formulation
Generic templates
ExaPF.AbstractVariable — TypeAbstractVariable
Variables corresponding to a particular formulation.
sourceExaPF.AbstractFormulation — TypeAbstractFormulation
Interface between the data and the mathemical formulation.
sourceExaPF.State — TypeState <: AbstractVariable
All variables $x$ depending on the variables Control $u$ through the non-linear equation $g(x, u) = 0$.
sourceExaPF.Control — TypeControl <: AbstractVariable
Independent variables $u$ used in the reduced-space formulation.
sourceStructure and variables
ExaPF.PolarForm — TypePolarForm{T, IT, VT, MT} <: AbstractPolarFormulation
Implement 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");
@@ -12,7 +12,7 @@
giving a mathematical formulation with:
#controls: 5
#states : 14
-
sourceExaPF.BlockPolarForm — TypeBlockPolarForm{T, IT, VT, MT} <: AbstractFormulation
Block polar formulation: duplicates k different polar models to evaluate them in parallel.
sourceExaPF.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")
+
sourceExaPF.BlockPolarForm — TypeBlockPolarForm{T, IT, VT, MT} <: AbstractFormulation
Block polar formulation: duplicates k different polar models to evaluate them in parallel.
sourceExaPF.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)
@@ -21,7 +21,7 @@
giving a mathematical formulation with:
#controls: 5
#states : 14
-
sourceExaPF.NetworkStack — TypeNetworkStack{VT,VD,MT} <: AbstractNetworkStack{VT}
+
sourceExaPF.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");
@@ -39,7 +39,7 @@
1.0
1.0
1.0
-
sourceExaPF.init! — Functioninit!(polar::PolarForm, stack::NetworkStack)
Set stack.input with the initial values specified in the base PS.PowerNetwork object.
sourceThe 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");
+
sourceExaPF.init! — Functioninit!(polar::PolarForm, stack::NetworkStack)
Set stack.input with the initial values specified in the base PS.PowerNetwork object.
sourceThe 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}:
@@ -48,7 +48,7 @@
3
20
21
-
sourcemapping(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");
+
sourcemapping(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}:
@@ -66,7 +66,7 @@
7
8
9
-
sourcePowerflow solver
ExaPF.run_pf — Functionrun_pf(
+
sourcePowerflow solver
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 network
Examples
julia> polar = ExaPF.load_polar("case9");
@@ -79,7 +79,7 @@
true
julia> conv.n_iterations
-4
sourceExaPF.nlsolve! — Functionnlsolve!(
+4
sourceExaPF.nlsolve! — Functionnlsolve!(
algo::NewtonRaphson,
jac::Jacobian,
stack::NetworkStack;
@@ -100,7 +100,7 @@
true
julia> conv.n_iterations
-4
sourceExaPF.NewtonRaphson — TypeNewtonRaphson <: AbstractNonLinearSolver
Newton-Raphson algorithm.
Attributes
maxiter::Int (default 20): maximum number of iterationstol::Float64 (default 1e-8): tolerance of the algorithmverbose::Int (default 0): verbosity level
sourceConstraints
The different parts of the polar formulation are implemented in the following AbstractExpression:
ExaPF.PolarBasis — TypePolarBasis{VI, MT} <: AbstractExpression
+4
sourceExaPF.NewtonRaphson — TypeNewtonRaphson <: AbstractNonLinearSolver
Newton-Raphson algorithm.
Attributes
maxiter::Int (default 20): maximum number of iterationstol::Float64 (default 1e-8): tolerance of the algorithmverbose::Int (default 0): verbosity level
sourceConstraints
The 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 \\
@@ -134,7 +134,7 @@
1.0
1.0
1.0
-
sourceExaPF.PowerFlowBalance — TypePowerFlowBalance{VT, MT}
+
sourceExaPF.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\} \\
@@ -167,7 +167,7 @@
julia> isapprox(powerflow(stack), zeros(14); atol=1e-8)
true
-
sourceExaPF.VoltageMagnitudeBounds — TypeVoltageMagnitudeBounds
Implement 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");
+
sourceExaPF.VoltageMagnitudeBounds — TypeVoltageMagnitudeBounds
Implement 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);
@@ -181,7 +181,7 @@
1.0
1.0
1.0
-
sourceExaPF.PowerGenerationBounds — TypePowerGenerationBounds{VT, MT}
+
sourceExaPF.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");
@@ -197,7 +197,7 @@
0.24069
0.144601
-0.03649
-
sourceExaPF.LineFlows — TypeLineFlows{VT, MT}
+
sourceExaPF.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);
@@ -226,7 +226,7 @@
2.67781
0.726668
0.215497
-
sourceObjective
The production costs is given in the AbstractExpression CostFunction:
ExaPF.CostFunction — TypeCostFunction{VT, MT} <: AutoDiff.AbstractExpression
+
sourceObjective
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);
@@ -236,5 +236,5 @@
julia> cost(stack)
1-element Vector{Float64}:
4509.0275
-
sourceComposition of expressions
The different expressions can be combined together in several different ways.
ExaPF.MultiExpressions — TypeMultiExpressions <: AbstractExpression
Implement 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)]
-
sourceExaPF.ComposedExpressions — TypeComposedExpressions{Expr1<:PolarBasis, Expr2} <: AbstractExpression
Implement 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.
sourceSettings
This document was generated with Documenter.jl version 0.27.25 on Saturday 2 September 2023. Using Julia version 1.8.5.
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 Thursday 21 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;
+Linear Solvers · ExaPF.jl Linear solvers
Description
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: RHS
Solve 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.
sourceDirect solvers
ExaPF wraps UMFPACK (shipped with SuiteSparse.jl) on the CPU, and CUSPARSE on CUDA device.
ExaPF.LinearSolvers.DirectSolver — TypeDirectSolver <: AbstractLinearSolver
Solve linear system $A x = y$ with direct linear algebra.
- On the CPU,
DirectSolver uses UMFPACK to solve the linear system - On CUDA GPU,
DirectSolver redirects the resolution to the method CUSOLVER.csrlsvqr
sourceIterative solvers
ExaPF.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$.
sourceExaPF.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$.
sourceExaPF.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$.
sourceExaPF.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$.
sourceExaPF.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.
sourceAvailable linear solvers can be queried with
ExaPF.LinearSolvers.list_solvers — Functionlist_solvers(::KernelAbstractions.Device)
List linear solvers available on current device.
sourcePreconditioning
To solve linear systems with iterative methods, ExaPF provides an implementation of a block-Jacobi preconditioner, portable on GPU.
ExaPF.LinearSolvers.AbstractPreconditioner — TypeAbstractPreconditioner
Preconditioners for the iterative solvers mostly focused on GPUs
sourceBlock-Jacobi preconditioner
ExaPF.LinearSolvers.BlockJacobiPreconditioner — TypeBlockJacobiPreconditioner
Overlapping-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 GPU
sourceExaPF.LinearSolvers.update — Functionfunction update(J::CuSparseMatrixCSR, p)
Update the preconditioner p from the sparse Jacobian J in CSR format for the GPU
- The dense blocks
cuJs are filled from the sparse Jacobian J - To a batch inversion of the dense blocks using CUBLAS
- Extract the preconditioner matrix
p.P from the dense blocks cuJs
sourcefunction 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.
sourceSettings
This document was generated with Documenter.jl version 0.27.25 on Saturday 2 September 2023. Using Julia version 1.8.5.
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 Thursday 21 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
+PowerSystem · ExaPF.jl PowerSystem
Description
ExaPF.PowerSystem.AbstractPowerSystem — TypeAbstractPowerSystem
First layer of the package. Store the topology of a given transmission network, including:
- the power injection at each bus ;
- the admittance matrix ;
- the default voltage at each bus.
Data are imported either from a matpower file, or a PSSE file.
sourceExaPF.PowerSystem.PowerNetwork — TypePowerNetwork <: AbstractPowerSystem
This 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.
sourceExaPF.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.
@@ -7,6 +7,6 @@
julia> net = PS.load_case("case1354_pegase", PS.PGLIB)
PowerNetwork object with:
Buses: 1354 (Slack: 1. PV: 259. PQ: 1094)
- Generators: 260.
sourceAPI Reference
Network elements
ExaPF.PowerSystem.AbstractNetworkElement — TypeAbstractNetworkElement
Abstraction for all physical elements being parts of a AbstractPowerSystem. Elements are divided in
- transmission lines (
Lines) - buses (
Buses) - generators (
Generators)
sourceList of elements:
ExaPF.PowerSystem.Buses — TypeBuses <: AbstractNetworkElement
Buses of a transmission network.
sourceExaPF.PowerSystem.Lines — TypeLines <: AbstractNetworkElement
Lines of a transmission network.
sourceExaPF.PowerSystem.Generators — TypeGenerators <: AbstractElement
Generators in a transmission network
sourceNetwork attributes
ExaPF.PowerSystem.AbstractNetworkAttribute — TypeAbstractNetworkAttribute
Attribute of a AbstractPowerSystem.
sourceList of attributes:
ExaPF.PowerSystem.NumberOfBuses — TypeNumberOfBuses <: AbstractNetworkAttribute
Number of buses in a AbstractPowerSystem.
sourceExaPF.PowerSystem.NumberOfLines — TypeNumberOfLines <: AbstractNetworkAttribute
Number of lines in a AbstractPowerSystem.
sourceExaPF.PowerSystem.NumberOfGenerators — TypeNumberOfGenerators <: AbstractNetworkAttribute
Number of generators in a AbstractPowerSystem.
sourceExaPF.PowerSystem.NumberOfPVBuses — TypeNumberOfPVBuses <: AbstractNetworkAttribute
Number of PV buses in a AbstractPowerSystem.
sourceExaPF.PowerSystem.NumberOfPQBuses — TypeNumberOfPQBuses <: AbstractNetworkAttribute
Number of PQ buses in a AbstractPowerSystem.
sourceExaPF.PowerSystem.NumberOfSlackBuses — TypeNumberOfSlackBuses <: AbstractNetworkAttribute
Number of slack buses in a AbstractPowerSystem.
sourceExaPF.PowerSystem.BaseMVA — TypeBaseMVA <: AbstractNetworkAttribute
Base MVA of the network.
sourceExaPF.PowerSystem.BusAdmittanceMatrix — TypeBusAdmittanceMatrix <: AbstractNetworkAttribute
Bus admittance matrix associated with the topology of the network.
sourceQuery the indexing of the different elements in a given network:
ExaPF.PowerSystem.PVIndexes — TypePVIndexes <: AbstractIndexing
Indexes of the PV buses in a AbstractPowerSystem.
sourceExaPF.PowerSystem.PQIndexes — TypePQIndexes <: AbstractIndexing
Indexes of the PQ buses in a AbstractPowerSystem.
sourceExaPF.PowerSystem.SlackIndexes — TypeSlackIndexes <: AbstractIndexing
Indexes of the slack buses in a AbstractPowerSystem.
sourceExaPF.PowerSystem.GeneratorIndexes — TypeGeneratorIndexes <: AbstractIndexing
Indexes of the generators in a AbstractPowerSystem.
sourceNetwork values
ExaPF.PowerSystem.AbstractNetworkValues — TypeAbstractNetworkValues
Numerical values attached to the different attributes of the network.
sourceList of values:
ExaPF.PowerSystem.VoltageMagnitude — TypeVoltageMagnitude <: AbstractNetworkValues
Magnitude |v| of the voltage v = |v| exp(i θ).
sourceExaPF.PowerSystem.VoltageAngle — TypeVoltageAngle <: AbstractNetworkValues
Angle θ of the voltage v = |v| exp(i θ).
sourceExaPF.PowerSystem.ActivePower — TypeActivePower <: AbstractNetworkValues
Active power P of the complex power S = P + iQ.
sourceExaPF.PowerSystem.ReactivePower — TypeReactivePower <: AbstractNetworkValues
Reactive power Q of the complex power S = P + iQ.
sourceExaPF.PowerSystem.ActiveLoad — TypeActiveLoad <: AbstractNetworkValues
Active load Pd at buses.
sourceExaPF.PowerSystem.ReactiveLoad — TypeReactiveLoad <: AbstractNetworkValues
Reactive load Qd at buses.
sourceFunction 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())
+ Generators: 260.
sourceAPI Reference
Network elements
ExaPF.PowerSystem.AbstractNetworkElement — TypeAbstractNetworkElement
Abstraction for all physical elements being parts of a AbstractPowerSystem. Elements are divided in
- transmission lines (
Lines) - buses (
Buses) - generators (
Generators)
sourceList of elements:
ExaPF.PowerSystem.Buses — TypeBuses <: AbstractNetworkElement
Buses of a transmission network.
sourceExaPF.PowerSystem.Lines — TypeLines <: AbstractNetworkElement
Lines of a transmission network.
sourceExaPF.PowerSystem.Generators — TypeGenerators <: AbstractElement
Generators in a transmission network
sourceNetwork attributes
ExaPF.PowerSystem.AbstractNetworkAttribute — TypeAbstractNetworkAttribute
Attribute of a AbstractPowerSystem.
sourceList of attributes:
ExaPF.PowerSystem.NumberOfBuses — TypeNumberOfBuses <: AbstractNetworkAttribute
Number of buses in a AbstractPowerSystem.
sourceExaPF.PowerSystem.NumberOfLines — TypeNumberOfLines <: AbstractNetworkAttribute
Number of lines in a AbstractPowerSystem.
sourceExaPF.PowerSystem.NumberOfGenerators — TypeNumberOfGenerators <: AbstractNetworkAttribute
Number of generators in a AbstractPowerSystem.
sourceExaPF.PowerSystem.NumberOfPVBuses — TypeNumberOfPVBuses <: AbstractNetworkAttribute
Number of PV buses in a AbstractPowerSystem.
sourceExaPF.PowerSystem.NumberOfPQBuses — TypeNumberOfPQBuses <: AbstractNetworkAttribute
Number of PQ buses in a AbstractPowerSystem.
sourceExaPF.PowerSystem.NumberOfSlackBuses — TypeNumberOfSlackBuses <: AbstractNetworkAttribute
Number of slack buses in a AbstractPowerSystem.
sourceExaPF.PowerSystem.BaseMVA — TypeBaseMVA <: AbstractNetworkAttribute
Base MVA of the network.
sourceExaPF.PowerSystem.BusAdmittanceMatrix — TypeBusAdmittanceMatrix <: AbstractNetworkAttribute
Bus admittance matrix associated with the topology of the network.
sourceQuery the indexing of the different elements in a given network:
ExaPF.PowerSystem.PVIndexes — TypePVIndexes <: AbstractIndexing
Indexes of the PV buses in a AbstractPowerSystem.
sourceExaPF.PowerSystem.PQIndexes — TypePQIndexes <: AbstractIndexing
Indexes of the PQ buses in a AbstractPowerSystem.
sourceExaPF.PowerSystem.SlackIndexes — TypeSlackIndexes <: AbstractIndexing
Indexes of the slack buses in a AbstractPowerSystem.
sourceExaPF.PowerSystem.GeneratorIndexes — TypeGeneratorIndexes <: AbstractIndexing
Indexes of the generators in a AbstractPowerSystem.
sourceNetwork values
ExaPF.PowerSystem.AbstractNetworkValues — TypeAbstractNetworkValues
Numerical values attached to the different attributes of the network.
sourceList of values:
ExaPF.PowerSystem.VoltageMagnitude — TypeVoltageMagnitude <: AbstractNetworkValues
Magnitude |v| of the voltage v = |v| exp(i θ).
sourceExaPF.PowerSystem.VoltageAngle — TypeVoltageAngle <: AbstractNetworkValues
Angle θ of the voltage v = |v| exp(i θ).
sourceExaPF.PowerSystem.ActivePower — TypeActivePower <: AbstractNetworkValues
Active power P of the complex power S = P + iQ.
sourceExaPF.PowerSystem.ReactivePower — TypeReactivePower <: AbstractNetworkValues
Reactive power Q of the complex power S = P + iQ.
sourceExaPF.PowerSystem.ActiveLoad — TypeActiveLoad <: AbstractNetworkValues
Active load Pd at buses.
sourceExaPF.PowerSystem.ReactiveLoad — TypeReactiveLoad <: AbstractNetworkValues
Reactive load Qd at buses.
sourceFunction 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())
-
sourceSettings
This document was generated with Documenter.jl version 0.27.25 on Saturday 2 September 2023. Using Julia version 1.8.5.
Settings
This document was generated with Documenter.jl version 1.0.1 on Thursday 21 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
+AutoDiff · ExaPF.jl AutoDiff
Overview
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
@@ -35,4 +35,4 @@
F[npq + i] += coef_cos * sin_val - coef_sin * cos_val
end
end
-end
These 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.
- 1Griewank, Andreas, and Andrea Walther. Evaluating derivatives: principles and techniques of algorithmic differentiation. Society for Industrial and Applied Mathematics, 2008.
Settings
This document was generated with Documenter.jl version 0.27.25 on Saturday 2 September 2023. Using Julia version 1.8.5.
These 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 Thursday 21 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)
+Benchmark · ExaPF.jl Benchmark
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.m
yields
KrylovBICGSTAB, CUDABackend, case300.m, 69.0, 3.57, 43.7, true
The first three fields are the settings of the benchmark run. They are followed by three timings in milliseconds:
- the time taken by the Newton-Raphson algorithm to solve the power flow,
- the timings for the Jacobian accumulation using AutoDiff,
- and the time for the linear solver, including the preconditioner.
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 0.27.25 on Saturday 2 September 2023. Using Julia version 1.8.5.
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 Thursday 21 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 \\ +
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 \\ @@ -173,4 +173,4 @@ 0.0 0.02340899999999987 0.007743999999999858 -
Settings
This document was generated with Documenter.jl version 0.27.25 on Saturday 2 September 2023. Using Julia version 1.8.5.
Settings
This document was generated with Documenter.jl version 1.0.1 on Thursday 21 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 0.27.25 on Saturday 2 September 2023. Using Julia version 1.8.5.
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 Thursday 21 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
+PowerSystem · ExaPF.jl PowerSystem
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
Description
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})) \,.
@@ -16,4 +16,4 @@
2
julia> PS.get(ps, PS.NumberOfSlackBuses())
-1
Settings
This document was generated with Documenter.jl version 0.27.25 on Saturday 2 September 2023. Using Julia version 1.8.5.
Settings
This document was generated with Documenter.jl version 1.0.1 on Thursday 21 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
+Quick start · ExaPF.jl Quick Start
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 KernelAbstractions
Then, we load ExaPF and its submodules with
using ExaPF
import ExaPF: AutoDiff
const PS = ExaPF.PowerSystem
@@ -15,12 +15,12 @@
#it 2: 5.88242e-01
#it 3: 4.88493e-03
#it 4: 1.39924e-06
-#it 5: 1.35671e-11
+#it 5: 1.52178e-11
Power flow has converged: true
* #iterations: 5
- * Time Jacobian (s) ........: 0.5833
- * Time linear solver (s) ...: 0.0210
- * Time total (s) ...........: 1.1181
Implicitly, ExaPF has just proceed to the following operations:
- instantiate automatically a starting point $x_0$ from the variables stored in
stack - instantiate the Jacobian of the powerflow equations using AutoDiff.
- solve the powerflow equations iteratively, using a Newton-Raphson algorithm.
This 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.
Detailed version
In what follows, we detail step by step the detailed procedure to solve the powerflow equations.
How to load a MATPOWER instance as a PowerNetwork object?
We start by importing a MATPOWER instance to a ExaPF.PowerSystem.PowerNetwork object:
julia> pf = PS.PowerNetwork(datafile)PowerNetwork object with:
+ * Time Jacobian (s) ........: 0.4023
+ * Time linear solver (s) ...: 0.0168
+ * Time total (s) ...........: 0.7970
Implicitly, ExaPF has just proceed to the following operations:
- instantiate automatically a starting point $x_0$ from the variables stored in
stack - instantiate the Jacobian of the powerflow equations using AutoDiff.
- solve the powerflow equations iteratively, using a Newton-Raphson algorithm.
This 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.
Detailed version
In what follows, we detail step by step the detailed procedure to solve the powerflow equations.
How to load a MATPOWER instance as a PowerNetwork object?
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())1354
julia> 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.
How to solve the powerflow equations?
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})) \,, \\
@@ -34,36 +34,36 @@
#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:
- an initial variable $x_0$
- a function to solve efficiently the linear system $(\nabla_x g_k) x_{k+1} = g(x_k, u)$
- a function to evaluate the Jacobian $\nabla_x g_k$
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
+⎡⢿⣷⣲⣷⣯⣷⣣⢯⣮⣌⢫⡿⣿⣲⢡⣾⡾⣶⣰⣷⡧⣾⣵⣿⢺⣵⣧⣙⣿⢿⣟⡮⢴⣗⣷⣶⣾⣿⡱⡧⎤
+⎢⢼⣾⣿⣿⣿⡺⣶⣿⣟⢷⣲⣿⣷⣿⣳⡟⢹⡿⣿⢿⢓⣻⡗⣾⢼⡿⡿⣆⣾⣾⣯⣗⢿⠃⡻⣿⡿⣟⢚⣧⎥
+⎢⢯⣿⣻⡻⢟⣵⢿⣟⣇⣾⣵⣿⡿⡷⣿⣿⢿⣿⢿⢟⣿⡽⢿⣾⢟⣾⣷⣫⣾⣿⣿⡾⣯⣾⣟⡿⡷⣿⣿⠿⎥
+⎢⡭⣞⣼⣿⣿⢷⡿⣯⣧⣿⢿⣾⣷⣷⡫⣧⣻⣻⣵⢾⡹⣻⡿⣿⢿⣿⣹⣾⣿⣿⣶⡟⣻⣌⣿⣮⣷⡏⣿⡇⎥
+⎢⡊⢿⢿⣝⣩⣽⣭⣿⣿⢟⣻⣿⡷⢿⢽⣯⣿⢾⣼⡻⢽⣼⢯⡻⣽⢿⡿⣟⣿⣿⠿⣿⣿⣽⡷⣿⣽⡯⣧⣿⎥
+⎢⣯⡶⣼⣾⣵⣿⣻⣷⣿⣾⡛⣬⣛⣿⢻⣽⠽⣯⣿⣳⣙⣿⣾⣝⣿⣴⡿⣿⢳⣽⣻⡯⣯⡿⢧⣻⣾⣏⣻⢗⎥
+⎢⢻⣻⣽⣿⢿⡯⢽⣿⣽⣏⣿⣼⡿⣯⣞⣷⡿⣹⡽⣗⣾⣟⢵⡭⣿⡯⣷⣹⣧⣻⢿⣷⣳⣳⠯⣿⣯⣷⣾⡝⎥
+⎢⣡⣶⣽⠾⣿⣿⠯⣮⡷⣷⣟⣶⢾⣽⣻⣾⣷⣷⢿⣿⣾⣳⣿⣿⢵⣾⢿⣾⣻⣦⢯⣞⣿⣿⣿⢿⣟⣿⣟⣜⎥
+⎢⢺⣯⣷⡶⣿⣷⣿⣺⣻⣟⡷⣧⣟⣫⢽⣿⡻⣮⣿⣳⢻⣿⣷⢿⢟⣞⣯⣷⣻⣾⣓⣩⣿⣟⢿⣾⣟⡮⣺⣿⎥
+⎢⢴⣾⣿⣟⣿⢗⣱⣟⣶⡻⢿⣻⢷⢯⣿⣷⢿⣻⢿⢗⣚⣞⣾⣒⣿⣷⣖⡷⣿⡿⡼⣿⣿⣾⣿⣯⣿⣖⢺⡹⎥
+⎢⣩⣯⣽⣰⣟⡿⣷⣪⣓⣷⣷⣼⣾⢿⢾⣻⣿⣶⣺⢼⡿⣯⣿⢿⣗⣝⣹⣿⢷⣷⡾⡷⢟⢿⣷⣲⡯⣞⣿⣿⎥
+⎢⣵⣿⣹⣭⣻⣷⣿⣯⣯⡳⣞⢿⡕⡷⣿⣿⣽⣟⢺⢻⣿⣟⣿⣿⢯⣟⣟⢶⣾⣻⣻⣿⣷⢯⣿⣟⡟⣿⡽⡻⎥
+⎢⢞⣶⣶⡷⣻⣵⣿⣷⣷⣟⢛⣿⡿⡿⣱⣷⣻⢵⢿⣿⣝⢽⣯⢷⣿⣿⣽⣟⣻⣿⢿⠏⣹⣦⡿⢶⣿⣏⡿⢇⎥
+⎢⣍⢻⠻⢯⡽⣻⣳⣾⣿⢯⣿⣯⣝⣻⣻⣷⢯⣿⢼⡽⣷⣾⢻⣝⣷⢿⣟⣽⡻⣭⣛⣽⢿⡿⡿⣳⣾⣗⣷⡻⎥
+⎢⣿⣟⣺⣿⣾⣿⣿⣿⣿⣿⣝⣶⣭⣻⠻⣾⣻⣾⣿⡿⢽⣷⣾⣻⣿⣾⡟⣮⣻⣾⣝⡷⢷⡟⣧⣾⣿⡯⢿⣲⎥
+⎢⡻⡽⢯⢿⣻⡿⣼⠿⣿⣧⡿⡾⢿⣷⣫⢷⡝⣸⣶⣯⢾⡯⣿⣾⡿⠗⣟⣼⢷⡽⣿⣿⣹⣿⣗⡿⣦⣷⣿⣎⎥
+⎢⢴⢷⠿⠓⣫⣿⡛⢾⣟⣿⣯⡿⢽⣺⣿⣿⣿⢿⣻⣿⣿⣕⡽⣟⠳⣾⣿⡷⣽⠷⣷⣾⣿⣿⣿⢿⣿⣿⣯⣶⎥
+⎢⢹⣿⣿⣮⣿⡽⡻⣿⣽⣯⣭⣳⣯⣧⣿⣟⣻⣷⡿⣿⢹⣻⣿⢿⢻⣏⢿⣫⣩⣿⣽⡽⣿⣟⣵⣿⣿⡷⣽⣷⎥
+⎢⣾⣿⣿⢯⣽⣯⡽⠿⡷⡿⡾⢿⢯⣿⣿⣽⡻⡽⢻⢿⣫⢯⣿⣭⡿⢿⢾⢿⡿⡿⢬⣿⣿⣿⢿⡿⣵⣿⣹⣝⎥
+⎣⠵⡮⠾⣴⣿⡟⠿⠿⣭⣿⢿⢞⣞⠿⣛⢽⣾⣾⣞⡲⣿⣿⣷⡫⠿⢏⣽⡻⢻⣳⡻⢿⢫⣿⢷⣿⣗⢾⣿⣿⎦
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.35671e-11
+#it 5: 1.52178e-11
Power flow has converged: true
* #iterations: 5
- * Time Jacobian (s) ........: 0.0040
- * Time linear solver (s) ...: 0.0193
- * Time total (s) ...........: 0.0236
Here, the algorithm solves the powerflow equations in 5 iterations. The algorithm modifies the values of stack inplace, to avoid any unnecessary memory allocations.
How to deport the computation on the GPU?
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))
+ * Time Jacobian (s) ........: 0.0045
+ * Time linear solver (s) ...: 0.0168
+ * Time total (s) ...........: 0.0216
Here, the algorithm solves the powerflow equations in 5 iterations. The algorithm modifies the values of stack inplace, to avoid any unnecessary memory allocations.
How to deport the computation on the GPU?
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
@@ -76,19 +76,19 @@
#it 2: 5.88242e-01
#it 3: 4.88493e-03
#it 4: 1.39924e-06
-#it 5: 1.23093e-11
+#it 5: 9.96622e-12
Power flow has converged: true
* #iterations: 5
- * Time Jacobian (s) ........: 3.1577
- * Time linear solver (s) ...: 0.3269
- * Time total (s) ...........: 4.5047
Note that we get exactly the same iterations as when we solve the power flow equations on the CPU.
How to solve the linear system with BICGSTAB?
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
+ * Time Jacobian (s) ........: 2.6187
+ * Time linear solver (s) ...: 0.3118
+ * Time total (s) ...........: 3.7690
Note that we get exactly the same iterations as when we solve the power flow equations on the CPU.
How to solve the linear system with BICGSTAB?
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: 1.35570e-11
+#it 5: 6.63821e-11
Power flow has converged: true
* #iterations: 5
- * Time Jacobian (s) ........: 0.0035
- * Time linear solver (s) ...: 6.8257
- * Time total (s) ...........: 6.8306
Settings
This document was generated with Documenter.jl version 0.27.25 on Saturday 2 September 2023. Using Julia version 1.8.5.
Settings
This document was generated with Documenter.jl version 1.0.1 on Thursday 21 September 2023. Using Julia version 1.9.3.
Settings
This document was generated with Documenter.jl version 0.27.25 on Saturday 2 September 2023. Using Julia version 1.8.5.
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());
+Power flow: batch evaluation · ExaPF.jl 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
@@ -21,26 +21,26 @@
#states : 14
Then, 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.12765746063391836
- 0.566960987817609
- 0.8524691406094105
- 0.12720379942550342
- 0.25327332442079364
- 0.5949761078937227
- 0.5558563498768199
- 0.26343616787518
- 0.38672875467598355
- 0.6053038499903975
+ 0.5941674359224426
+ 0.8522137692943615
+ 0.898905599393412
+ 0.7706263989403442
+ 0.10242552305040731
+ 0.1093694659440686
+ 0.24382786436620574
+ 0.2752153853589945
+ 0.38808902745269813
+ 0.010479031681545958
⋮
- 0.6609679827337362
- 0.5498645637015432
- 0.270136908048468
- 0.1434018453437842
- 0.2560352432323393
- 0.04200672968921837
- 0.4201631619773307
- 0.7724781105966682
- 0.6213140678547853
The 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}:
+ 0.43651484710096244
+ 0.33871529924162136
+ 0.04994653600628185
+ 0.004424390430899927
+ 0.2871564676452214
+ 0.14405974048525894
+ 0.8920382389401823
+ 0.47054505640394984
+ 0.4780947294409694
The 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
@@ -50,92 +50,92 @@
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
Only 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.623778 0.546067 0.706786 0.0100577 … 0.11351 0.472291 0.580118
- 0.169438 0.603691 0.122567 0.511803 0.724155 0.297848 0.774841
- 0.469862 0.827131 0.858884 0.758878 0.16326 0.754959 0.384554
- 0.0335776 0.832894 0.326203 0.786805 0.12509 0.964754 0.0576321
- 0.57805 0.953554 0.682215 0.612529 0.908971 0.360122 0.668147
- 0.417127 0.89729 0.186815 0.249881 … 0.895317 0.862494 0.230494
- 0.641264 0.369938 0.0298789 0.856954 0.115591 0.83833 0.340071
- 0.0692794 0.646089 0.406713 0.534395 0.179417 0.524836 0.694479
- 0.576422 0.202107 0.806112 0.107095 0.715618 0.108914 0.827546
Evaluate 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);
+ 0.81872 0.58149 0.993387 0.296823 … 0.985247 0.390188 0.233188
+ 0.176536 0.79199 0.362258 0.0995252 0.871698 0.373152 0.295096
+ 0.705494 0.891877 0.942154 0.562401 0.807404 0.499716 0.348183
+ 0.749389 0.368638 0.161537 0.209402 0.112298 0.452152 0.501078
+ 0.292929 0.378692 0.345064 0.829092 0.412382 0.198977 0.616138
+ 0.504299 0.95563 0.736151 0.793169 … 0.765123 0.163487 0.995522
+ 0.616872 0.836865 0.494355 0.370883 0.0303673 0.986578 0.878993
+ 0.96446 0.881766 0.488783 0.295786 0.505388 0.94591 0.454449
+ 0.629922 0.426857 0.660948 0.0796323 0.631545 0.637177 0.269932
Evaluate 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.46056 -1.02631 … -0.905845 -1.33215 -0.855159
- -0.380138 -0.0228686 -0.68674 -0.0950409 -0.465446
- 0.0335776 0.832894 0.12509 0.964754 0.0576321
- 0.57805 0.953554 0.908971 0.360122 0.668147
- 0.417127 0.89729 0.895317 0.862494 0.230494
- 0.641264 0.369938 … 0.115591 0.83833 0.340071
- 0.0692794 0.646089 0.179417 0.524836 0.694479
- 0.576422 0.202107 0.715618 0.108914 0.827546
- -0.0397962 0.00187016 0.574151 0.651935 -0.0235982
- -0.00472668 0.254469 -0.207384 0.210021 -0.00196476
- 0.311476 0.0389193 … 0.539286 0.0769768 -0.241493
- 0.376856 0.449859 -0.0912412 0.778293 0.241163
- 0.0359362 -0.105802 -0.200044 0.63735 0.544978
- 0.145729 -0.167072 0.451261 -0.144298 0.380314
We 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());
+ -1.45346 -0.83801 -1.26774 … -0.758302 -1.25685 -1.3349
+ -0.144506 0.0418768 0.0921542 -0.0425958 -0.350284 -0.501817
+ 0.749389 0.368638 0.161537 0.112298 0.452152 0.501078
+ 0.292929 0.378692 0.345064 0.412382 0.198977 0.616138
+ 0.504299 0.95563 0.736151 0.765123 0.163487 0.995522
+ 0.616872 0.836865 0.494355 … 0.0303673 0.986578 0.878993
+ 0.96446 0.881766 0.488783 0.505388 0.94591 0.454449
+ 0.629922 0.426857 0.660948 0.631545 0.637177 0.269932
+ 0.603626 0.107701 0.620703 0.400412 0.0744612 -0.162576
+ -0.155574 -0.254843 0.498071 0.397487 0.361444 0.0291565
+ -0.174131 0.715848 0.36658 … 0.294529 0.543651 -0.13944
+ 0.0648279 -0.164528 0.0793475 0.0566474 -0.159946 0.713038
+ 0.0477154 -0.0382624 0.647422 0.214784 0.115698 0.243045
+ 0.147089 0.123457 0.352825 0.638275 -0.00369888 0.237095
We 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.J
140×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);
+⎡⡱⣮⡲⣞⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎤
+⎢⣸⢮⣻⣾⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
+⎢⠀⠀⠀⠀⡱⣮⡲⣞⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
+⎢⠀⠀⠀⠀⣸⢮⣻⣾⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
+⎢⠀⠀⠀⠀⠀⠀⠀⠀⡵⣯⡶⣝⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
+⎢⠀⠀⠀⠀⠀⠀⠀⠀⣜⢯⡻⣮⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
+⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡵⣯⡶⣝⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
+⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣜⢯⡻⣮⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
+⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡵⣯⡶⣝⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
+⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣜⢯⡻⣮⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
+⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡵⣯⣺⣜⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
+⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣚⢾⡻⣮⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
+⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡵⣯⣺⣜⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
+⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣚⢾⡻⣮⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
+⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡵⣯⣺⣜⠀⠀⠀⠀⠀⠀⠀⠀⎥
+⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣚⢾⡻⣮⠀⠀⠀⠀⠀⠀⠀⠀⎥
+⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡵⣯⡺⣕⠀⠀⠀⠀⎥
+⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢞⢮⡿⣯⠀⠀⠀⠀⎥
+⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡵⣯⡺⣕⎥
+⎣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢞⢮⡿⣯⎦
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!(
+⎡⡱⣮⡲⣞⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎤
+⎢⣸⢮⣻⣾⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
+⎢⠀⠀⠀⠀⡱⣮⡲⣞⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
+⎢⠀⠀⠀⠀⣸⢮⣻⣾⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
+⎢⠀⠀⠀⠀⠀⠀⠀⠀⡵⣯⡶⣝⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
+⎢⠀⠀⠀⠀⠀⠀⠀⠀⣜⢯⡻⣮⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
+⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡵⣯⡶⣝⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
+⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣜⢯⡻⣮⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
+⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡵⣯⡶⣝⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
+⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣜⢯⡻⣮⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
+⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡵⣯⣺⣜⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
+⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣚⢾⡻⣮⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
+⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡵⣯⣺⣜⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
+⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣚⢾⡻⣮⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
+⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡵⣯⣺⣜⠀⠀⠀⠀⠀⠀⠀⠀⎥
+⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣚⢾⡻⣮⠀⠀⠀⠀⠀⠀⠀⠀⎥
+⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡵⣯⡺⣕⠀⠀⠀⠀⎥
+⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢞⢮⡿⣯⠀⠀⠀⠀⎥
+⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡵⣯⡺⣕⎥
+⎣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢞⢮⡿⣯⎦
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
+ * #iterations: 4
* Time Jacobian (s) ........: 0.0002
- * Time linear solver (s) ...: 0.0006
- * Time total (s) ...........: 0.0010
+ * Time linear solver (s) ...: 0.0004
+ * Time total (s) ...........: 0.0007
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.984247 0.978034 0.988546 0.956992 0.956072 0.943995 0.977201
- 0.970967 0.937706 0.9936 0.964167 0.959314 0.918957 0.972038
- 0.970115 0.965811 0.99426 0.978848 … 0.969263 0.941392 0.990103
- 0.951801 0.949511 0.963229 0.974353 0.983533 0.886098 0.957444
- 0.976508 0.979187 0.985071 0.973009 0.98936 0.925285 0.956566
- 0.96449 0.986155 0.985193 0.939804 0.931591 0.944245 0.938306
Solve 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
+ 0.966309 0.972902 0.921606 0.964459 0.932831 0.975283 0.984083
+ 0.979567 0.964036 0.88175 0.957589 0.902597 0.943513 0.966502
+ 0.992338 0.948066 0.937069 0.957202 … 0.952171 0.962881 0.969966
+ 0.979551 0.964292 0.925649 0.924367 0.949875 0.973275 0.920203
+ 0.9821 0.973861 0.933158 0.932569 0.955048 0.97987 0.954013
+ 0.956651 0.956936 0.894466 0.901908 0.891719 0.970809 0.957078
Solve 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);
@@ -152,8 +152,8 @@
blk_stack_gpu;
linear_solver=rf_solver,
)
Power flow has converged: true
- * #iterations: 5
- * Time Jacobian (s) ........: 0.0016
- * Time linear solver (s) ...: 0.0636
- * Time total (s) ...........: 0.0666
-
Settings
This document was generated with Documenter.jl version 0.27.25 on Saturday 2 September 2023. Using Julia version 1.8.5.
Settings
This document was generated with Documenter.jl version 1.0.1 on Thursday 21 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
+Power flow: direct solver · ExaPF.jl Direct solvers for power flow
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.
UMFPACK (CPU, default)
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))
@@ -15,16 +15,16 @@
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.3883
- * Time linear solver (s) ...: 0.4640
- * Time total (s) ...........: 1.8579
+ * Time Jacobian (s) ........: 1.1515
+ * Time linear solver (s) ...: 0.3227
+ * Time total (s) ...........: 1.4793
KLU (CPU)
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} @0x00000000063589c0, Ptr{Nothing} @0x000000004f648f10, 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.25125875502204, 88.40062041224645, 995.2128749049234, -479.4136240982223, -138.29482949968116, -145.959631463846, -230.91975549270293, 22.972728768573436 … -480.47843944261814, -39.711370520685634, 520.1642960016138, -3360.900457400195, -208.11090119997567, 3568.982258669562, -39.404035452166724, 39.405224228971264, -242.8021538639655, 242.8019609359876]))
and pass it to nlsolve!:
ExaPF.init!(polar, stack) # reinit stack
+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} @0x0000000051ef7890, Ptr{Nothing} @0x000000004b89b190, 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.2602
- * Time linear solver (s) ...: 0.1734
- * Time total (s) ...........: 0.4387
+ * Time Jacobian (s) ........: 0.2700
+ * Time linear solver (s) ...: 0.0335
+ * Time total (s) ...........: 0.3084
We observe KLU reduces considerably the time spent in the linear solver.
cusolverRF (CUDA)
cusolverRF is an efficient LU refactorization routine implemented in CUDA. It is wrapped in Julia inside the package CUSOLVERRF.jl:
using CUSOLVERRF
The 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)
@@ -34,7 +34,7 @@
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.3278
- * Time linear solver (s) ...: 0.0828
- * Time total (s) ...........: 1.4121
-
Settings
This document was generated with Documenter.jl version 0.27.25 on Saturday 2 September 2023. Using Julia version 1.8.5.
Settings
This document was generated with Documenter.jl version 1.0.1 on Thursday 21 September 2023. Using Julia version 1.9.3.