From dc7ecdbb709807cf7781781dd16c34406e1f97a1 Mon Sep 17 00:00:00 2001 From: "Documenter.jl" Date: Tue, 5 Dec 2023 12:39:11 -0600 Subject: [PATCH] build based on 165d218 --- previews/PR284/.documenter-siteinfo.json | 2 +- previews/PR284/artifact/index.html | 2 +- previews/PR284/index.html | 2 +- previews/PR284/lib/autodiff/index.html | 10 +- previews/PR284/lib/formulations/index.html | 30 ++--- previews/PR284/lib/linearsolver/index.html | 18 +-- previews/PR284/lib/powersystem/index.html | 6 +- previews/PR284/man/autodiff/index.html | 2 +- previews/PR284/man/benchmark/index.html | 4 +- previews/PR284/man/formulations/index.html | 2 +- previews/PR284/man/linearsolver/index.html | 2 +- previews/PR284/man/powersystem/index.html | 2 +- previews/PR284/quickstart/index.html | 38 +++--- previews/PR284/search_index.js | 2 +- .../tutorials/batch_evaluation/index.html | 112 +++++++++--------- .../PR284/tutorials/direct_solver/index.html | 34 +++--- 16 files changed, 131 insertions(+), 137 deletions(-) diff --git a/previews/PR284/.documenter-siteinfo.json b/previews/PR284/.documenter-siteinfo.json index 11580ac0..f2dda870 100644 --- a/previews/PR284/.documenter-siteinfo.json +++ b/previews/PR284/.documenter-siteinfo.json @@ -1 +1 @@ -{"documenter":{"julia_version":"1.9.4","generation_timestamp":"2023-12-04T10:25:46","documenter_version":"1.2.1"}} \ No newline at end of file +{"documenter":{"julia_version":"1.9.4","generation_timestamp":"2023-12-05T12:39:05","documenter_version":"1.2.1"}} \ No newline at end of file diff --git a/previews/PR284/artifact/index.html b/previews/PR284/artifact/index.html index ab27ceae..fc2cb31a 100644 --- a/previews/PR284/artifact/index.html +++ b/previews/PR284/artifact/index.html @@ -1,2 +1,2 @@ -ExaData Artifact · ExaPF.jl

ExaData Artifact

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.

+ExaData Artifact · ExaPF.jl

ExaData Artifact

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.

diff --git a/previews/PR284/index.html b/previews/PR284/index.html index 9cb39df1..31af7685 100644 --- a/previews/PR284/index.html +++ b/previews/PR284/index.html @@ -1,2 +1,2 @@ -Home · ExaPF.jl

ExaPF

ExaPF.jl is a package to solve the power flow problem on upcoming exascale architectures. The code has been designed to be:

  1. Portable: Targeting exascale architectures implies a focus on graphics processing units (GPUs) as these systems lack substantial computational performance through classical CPUs.
  2. Differentiable: All the expressions implemented in ExaPF are fully compatible with ForwardDiff.jl, and routines are provided to extract first- and second-order derivatives to solve efficiently power flow and optimal power flow problems.

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 provide extensions that leverage the packages CUDA.jl, [AMDGPU.jl]((https://github.com/JuliaGPU/AMDGPU.jl), and KernelAbstractions.jl to make ExaPF portable across GPU architectures.

Table of contents

Manual

Library

Artifact

Funding

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.

+Home · ExaPF.jl

ExaPF

ExaPF.jl is a package to solve the power flow problem on upcoming exascale architectures. The code has been designed to be:

  1. Portable: Targeting exascale architectures implies a focus on graphics processing units (GPUs) as these systems lack substantial computational performance through classical CPUs.
  2. Differentiable: All the expressions implemented in ExaPF are fully compatible with ForwardDiff.jl, and routines are provided to extract first- and second-order derivatives to solve efficiently power flow and optimal power flow problems.

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 provide extensions that leverage the packages CUDA.jl, [AMDGPU.jl]((https://github.com/JuliaGPU/AMDGPU.jl), and KernelAbstractions.jl to make ExaPF portable across GPU architectures.

Table of contents

Manual

Library

Artifact

Funding

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.

diff --git a/previews/PR284/lib/autodiff/index.html b/previews/PR284/lib/autodiff/index.html index 385febaa..404f1aec 100644 --- a/previews/PR284/lib/autodiff/index.html +++ b/previews/PR284/lib/autodiff/index.html @@ -1,13 +1,13 @@ -AutoDiff · ExaPF.jl

AutoDiff

Variables

Expressions

ExaPF.AutoDiff.AbstractExpression — Type
AbstractExpression

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
+AutoDiff · ExaPF.jl
AutoDiff
Variables
Expressions
ExaPF.AutoDiff.AbstractExpression — Type
AbstractExpression
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.
source
ExaPF.AutoDiff.adjoint! — Function
adjoint!(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!.
source
First and second-order derivatives
Utils
ExaPF.AutoDiff.seed! — Function
seed!(
+
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.
source
ExaPF.AutoDiff.adjoint! — Function
adjoint!(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!.
source
First and second-order derivatives
Utils
ExaPF.AutoDiff.seed! — Function
seed!(
     H::AbstractHessianProd,
     v::AbstractVector{T},
-) where {T}
Seed the duals with v to compute the Hessian vector product $λ^⊤ H v$.
source
ExaPF.AutoDiff.seed_coloring! — Function
seed_coloring!(
     M::Union{AbstractJacobian, AbstractFullHessian}
     coloring::AbstractVector,
-)
Seed the duals with the coloring based seeds to compute the Jacobian or Hessian $M$.
source
ExaPF.AutoDiff.partials! — Function
partials!(jac::AbstractJacobian)
Extract partials from Jacobian jac in jac.J.
source
partials!(hess::AbstractFullHessian)
Extract partials from Hessian hess into hess.H.
source
ExaPF.AutoDiff.partials! — Function
partials!(jac::AbstractJacobian)
Extract partials from Jacobian jac in jac.J.
source
partials!(hess::AbstractFullHessian)
Extract partials from Hessian hess into hess.H.
source
ExaPF.AutoDiff.set_value! — Function
set_value!(
     jac,
     primals::AbstractVector{T}
-) where {T}
Set values of ForwardDiff.Dual numbers in jac to primals.
source

+) where {T}

Set values of ForwardDiff.Dual numbers in jac to primals.

source
diff --git a/previews/PR284/lib/formulations/index.html b/previews/PR284/lib/formulations/index.html index 4fc52664..30ef9c67 100644 --- a/previews/PR284/lib/formulations/index.html +++ b/previews/PR284/lib/formulations/index.html @@ -1,5 +1,5 @@ -Polar formulation · ExaPF.jl

Polar formulation

Generic templates

ExaPF.State — Type
State <: AbstractVariable

All variables $x$ depending on the variables Control $u$ through the non-linear equation $g(x, u) = 0$.

source
ExaPF.Control — Type
Control <: AbstractVariable

Independent variables $u$ used in the reduced-space formulation.

source

Structure and variables

ExaPF.PolarForm — Type
PolarForm{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;
+Polar formulation · ExaPF.jl
Polar formulation
Generic templates
ExaPF.State — Type
State <: AbstractVariable
All variables $x$ depending on the variables Control $u$ through the non-linear equation $g(x, u) = 0$.
source
ExaPF.Control — Type
Control <: AbstractVariable
Independent variables $u$ used in the reduced-space formulation.
source
Structure and variables
ExaPF.PolarForm — Type
PolarForm{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
-
source
ExaPF.BlockPolarForm — Type
BlockPolarForm{T, IT, VT, MT} <: AbstractFormulation
Block polar formulation: duplicates k different polar models to evaluate them in parallel.
source
ExaPF.load_polar — Function
load_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")
+
source
ExaPF.BlockPolarForm — Type
BlockPolarForm{T, IT, VT, MT} <: AbstractFormulation
Block polar formulation: duplicates k different polar models to evaluate them in parallel.
source
ExaPF.load_polar — Function
load_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
-
source
ExaPF.NetworkStack — Type
NetworkStack{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
-
source
The state and the control are defined as mapping:
ExaPF.mapping — Function
mapping(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");
+
source
The state and the control are defined as mapping:
ExaPF.mapping — Function
mapping(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
-
source
mapping(polar::PolarForm, ::State)
Return the mapping associated to the State() in NetworkStack according to the polar formulation PolarForm.
Examples
julia> polar = ExaPF.load_polar("case9");
+
source
mapping(polar::PolarForm, ::State)
Return the mapping associated to the State() in NetworkStack according to the polar formulation PolarForm.
Examples
julia> polar = ExaPF.load_polar("case9");
 
 julia> mapu = ExaPF.mapping(polar, State())
 14-element Vector{Int64}:
@@ -66,7 +66,7 @@
   7
   8
   9
-
source
Powerflow solver
Powerflow solver
ExaPF.run_pf — Function
run_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 equation
  • stack::NetworkStack: initial values in the network
Examples
julia> polar = ExaPF.load_polar("case9");
@@ -79,7 +79,7 @@
 true
 
 julia> conv.n_iterations
-4
source
ExaPF.nlsolve! — Function
nlsolve!(
     algo::NewtonRaphson,
     jac::Jacobian,
     stack::NetworkStack;
@@ -100,7 +100,7 @@
 true
 
 julia> conv.n_iterations
-4
source
ExaPF.NewtonRaphson — Type
NewtonRaphson <: AbstractNonLinearSolver
Newton-Raphson algorithm.
Attributes
  • maxiter::Int (default 20): maximum number of iterations
  • tol::Float64 (default 1e-8): tolerance of the algorithm
  • verbose::Int (default 0): verbosity level
source
Constraints
The different parts of the polar formulation are implemented in the following AbstractExpression:
ExaPF.NewtonRaphson — Type
NewtonRaphson <: AbstractNonLinearSolver
Newton-Raphson algorithm.
Attributes
  • maxiter::Int (default 20): maximum number of iterations
  • tol::Float64 (default 1e-8): tolerance of the algorithm
  • verbose::Int (default 0): verbosity level
source
Constraints
The different parts of the polar formulation are implemented in the following AbstractExpression:
ExaPF.PolarBasis — Type
PolarBasis{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
-
source
ExaPF.PowerFlowBalance — Type
PowerFlowBalance{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 -

source
ExaPF.VoltageMagnitudeBounds — Type
VoltageMagnitudeBounds

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");
+
source
ExaPF.VoltageMagnitudeBounds — Type
VoltageMagnitudeBounds

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
-
source
ExaPF.PowerGenerationBounds — Type
PowerGenerationBounds{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
-
source
ExaPF.LineFlows — Type
LineFlows{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
-
source

Objective

The production costs is given in the AbstractExpression CostFunction:

Objective

The production costs is given in the AbstractExpression CostFunction:

ExaPF.CostFunction — Type
CostFunction{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
-
source

Composition of expressions

The different expressions can be combined together in several different ways.

ExaPF.MultiExpressions — Type
MultiExpressions <: 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)]
-
source
ExaPF.ComposedExpressions — Type
ComposedExpressions{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.

source
+source

Composition of expressions

The different expressions can be combined together in several different ways.

ExaPF.MultiExpressions — Type
MultiExpressions <: 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)]
+
source
ExaPF.ComposedExpressions — Type
ComposedExpressions{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.

source
diff --git a/previews/PR284/lib/linearsolver/index.html b/previews/PR284/lib/linearsolver/index.html index 4baabbcb..53302202 100644 --- a/previews/PR284/lib/linearsolver/index.html +++ b/previews/PR284/lib/linearsolver/index.html @@ -1,16 +1,10 @@ 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!.

ExaPF.LinearSolvers.ldiv! — Function
ldiv!(solver, x, A, y)
-ldiv!(solver, x, y)
  • solver::AbstractLinearSolver: linear solver to solve the system
  • x::AbstractVector: Solution
  • A::AbstractMatrix: Input matrix
  • y::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.

source

Direct solvers

ExaPF wraps UMFPACK (shipped with SuiteSparse.jl) on the CPU, and CUSPARSE on CUDA device.

ExaPF.LinearSolvers.DirectSolver — Type
DirectSolver <: 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
source

Iterative solvers

ExaPF.LinearSolvers.KrylovBICGSTAB — Type
KrylovBICGSTAB <: 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$.

source
ExaPF.LinearSolvers.DQGMRES — Type
DQGMRES <: AbstractIterativeLinearSolver
-DQGMRES(precond; verbose=false, memory=4)

Wrap Krylov.jl's DQGMRES algorithm to solve iteratively the linear system $A x = y$.

source
ExaPF.LinearSolvers.BICGSTAB — Type
BICGSTAB <: AbstractIterativeLinearSolver
-BICGSTAB(precond; maxiter=2_000, tol=1e-8, verbose=false)

Custom BICGSTAB implementation to solve iteratively the linear system $A x = y$.

source
ExaPF.LinearSolvers.EigenBICGSTAB — Type
EigenBICGSTAB <: 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$.

source

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 — Function
bicgstab(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.

source

Available linear solvers can be queried with

A default solver is provided for each vendor backend.

Block-Krylov solvers

ExaPF.jl provides an implementation of block-GMRES to solve linear systems with multiple righ-hand sides.

ExaPF.LinearSolvers.BlockGmresSolver — Type

Type for storing the vectors required by the in-place version of BLOCK-GMRES.

The outer constructors

solver = BlockGmresSolver(m, n, p, memory, SV, SM)
-solver = BlockGmresSolver(A, B; memory=5)

may be used in order to create these vectors. memory is set to div(n,p) if the value given is larger than div(n,p).

source
ExaPF.LinearSolvers.block_gmres — Function
(X, stats) = block_gmres(A, B; X0::AbstractMatrix{FC}; memory::Int=20, M=I, N=I,
+ldiv!(solver, x, y)
  • solver::AbstractLinearSolver: linear solver to solve the system
  • x::AbstractVector: Solution
  • A::AbstractMatrix: Input matrix
  • y::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.

source

Direct solvers

ExaPF wraps UMFPACK (shipped with SuiteSparse.jl) on the CPU, and CUSPARSE on CUDA device.

ExaPF.LinearSolvers.DirectSolver — Type
DirectSolver <: 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
source

Iterative solvers

ExaPF.LinearSolvers.Bicgstab — Type
Bicgstab <: AbstractIterativeLinearSolver
+Bicgstab(precond; verbose=0, rtol=1e-10, atol=1e-10)

Wrap Krylov.jl's BICGSTAB algorithm to solve iteratively the linear system $A x = y$.

source
ExaPF.LinearSolvers.Dqgmres — Type
Dqgmres <: AbstractIterativeLinearSolver
+Dqgmres(precond; verbose=false, memory=4)

Wrap Krylov.jl's Dqgmres algorithm to solve iteratively the linear system $A x = y$.

source

Available linear solvers can be queried with

A default solver is provided for each vendor backend.

Block-Krylov solvers

ExaPF.jl provides an implementation of block-GMRES to solve linear systems with multiple righ-hand sides.

ExaPF.LinearSolvers.BlockGmresSolver — Type

Type for storing the vectors required by the in-place version of BLOCK-GMRES.

The outer constructors

solver = BlockGmresSolver(m, n, p, memory, SV, SM)
+solver = BlockGmresSolver(A, B; memory=5)

may be used in order to create these vectors. memory is set to div(n,p) if the value given is larger than div(n,p).

source
ExaPF.LinearSolvers.block_gmres — Function
(X, stats) = block_gmres(A, B; X0::AbstractMatrix{FC}; memory::Int=20, M=I, N=I,
                          ldiv::Bool=false, restart::Bool=false, reorthogonalization::Bool=false,
                          atol::T = √eps(T), rtol::T=√eps(T), itmax::Int=0,
-                         timemax::Float64=Inf, verbose::Int=0, history::Bool=false)

T is an AbstractFloat such as Float32, Float64 or BigFloat. FC is T or Complex{T}.

(X, stats) = block_gmres(A, B, X0::AbstractVector; kwargs...)

GMRES can be warm-started from an initial guess X0 where kwargs are the same keyword arguments as above.

Solve the linear system AX = B of size n with p right-hand sides using block-GMRES.

Input arguments

  • A: a linear operator that models a matrix of dimension n;
  • B: a matrix of size n × p.

Optional argument

  • X0: a matrix of size n × p that represents an initial guess of the solution X.

Keyword arguments

  • memory: if restart = true, the restarted version block-GMRES(k) is used with k = memory. If restart = false, the parameter memory should be used as a hint of the number of iterations to limit dynamic memory allocations. Additional storage will be allocated if the number of iterations exceeds memory;
  • M: linear operator that models a nonsingular matrix of size n used for left preconditioning;
  • N: linear operator that models a nonsingular matrix of size n used for right preconditioning;
  • ldiv: define whether the preconditioners use ldiv! or mul!;
  • restart: restart the method after memory iterations;
  • reorthogonalization: reorthogonalize the new matrices of the block-Krylov basis against all previous matrix;
  • atol: absolute stopping tolerance based on the residual norm;
  • rtol: relative stopping tolerance based on the residual norm;
  • itmax: the maximum number of iterations. If itmax=0, the default number of iterations is set to 2 * div(n,p);
  • timemax: the time limit in seconds;
  • verbose: additional details can be displayed if verbose mode is enabled (verbose > 0). Information will be displayed every verbose iterations;
  • history: collect additional statistics on the run such as residual norms.

Output arguments

  • x: a dense matrix of size n × p;
  • stats: statistics collected on the run in a BlockGmresStats.
source
+ timemax::Float64=Inf, verbose::Int=0, history::Bool=false)

T is an AbstractFloat such as Float32, Float64 or BigFloat. FC is T or Complex{T}.

(X, stats) = block_gmres(A, B, X0::AbstractVector; kwargs...)

GMRES can be warm-started from an initial guess X0 where kwargs are the same keyword arguments as above.

Solve the linear system AX = B of size n with p right-hand sides using block-GMRES.

Input arguments

Optional argument

Keyword arguments

Output arguments

source
ExaPF.LinearSolvers.block_gmres! — Function
solver = block_gmres!(solver::BlockGmresSolver, B; kwargs...)
+solver = block_gmres!(solver::BlockGmresSolver, B, X0; kwargs...)

where kwargs are keyword arguments of block_gmres.

See BlockGmresSolver for more details about the solver.

source
diff --git a/previews/PR284/lib/powersystem/index.html b/previews/PR284/lib/powersystem/index.html index ef688e90..2b158180 100644 --- a/previews/PR284/lib/powersystem/index.html +++ b/previews/PR284/lib/powersystem/index.html @@ -1,5 +1,5 @@ -PowerSystem · ExaPF.jl

PowerSystem

Description

ExaPF.PowerSystem.AbstractPowerSystem — Type
AbstractPowerSystem

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.

source
ExaPF.PowerSystem.PowerNetwork — Type
PowerNetwork <: 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.

source
ExaPF.PowerSystem.load_case — Function
load_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 — Type
AbstractPowerSystem
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.
source
ExaPF.PowerSystem.PowerNetwork — Type
PowerNetwork <: 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.
source
ExaPF.PowerSystem.load_case — Function
load_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.
source
API Reference
Network elements
ExaPF.PowerSystem.AbstractNetworkElement — Type
AbstractNetworkElement
Abstraction for all physical elements being parts of a AbstractPowerSystem. Elements are divided in
  • transmission lines (Lines)
  • buses (Buses)
  • generators (Generators)
source
List of elements:
Network attributes
List of attributes:
Query the indexing of the different elements in a given network:
Network values
List of values:
Function to get the range of a given value:
ExaPF.PowerSystem.bounds — Function
bounds(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.
source
API Reference
Network elements
ExaPF.PowerSystem.AbstractNetworkElement — Type
AbstractNetworkElement
Abstraction for all physical elements being parts of a AbstractPowerSystem. Elements are divided in
  • transmission lines (Lines)
  • buses (Buses)
  • generators (Generators)
source
List of elements:
Network attributes
List of attributes:
Query the indexing of the different elements in a given network:
Network values
List of values:
Function to get the range of a given value:
ExaPF.PowerSystem.bounds — Function
bounds(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())
-
source

+
source
diff --git a/previews/PR284/man/autodiff/index.html b/previews/PR284/man/autodiff/index.html index 04b8d6e7..81eded40 100644 --- a/previews/PR284/man/autodiff/index.html +++ b/previews/PR284/man/autodiff/index.html @@ -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

SIMD AD for point-wise vector product \\label{fig:simd}

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.

+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

SIMD AD for point-wise vector product \\label{fig:simd}

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.

diff --git a/previews/PR284/man/benchmark/index.html b/previews/PR284/man/benchmark/index.html index 65359035..31b9d30d 100644 --- a/previews/PR284/man/benchmark/index.html +++ b/previews/PR284/man/benchmark/index.html @@ -1,5 +1,5 @@ -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:

  1. the time taken by the Newton-Raphson algorithm to solve the power flow,
  2. the timings for the Jacobian accumulation using AutoDiff,
  3. 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)
+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. Bicgstab), 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 Bicgstab CUDABackend case300.m
yields
Bicgstab, 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:
  1. the time taken by the Newton-Raphson algorithm to solve the power flow,
  2. the timings for the Jacobian accumulation using AutoDiff,
  3. 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.

+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.

diff --git a/previews/PR284/man/formulations/index.html b/previews/PR284/man/formulations/index.html index 9dede310..b1218b4f 100644 --- a/previews/PR284/man/formulations/index.html +++ b/previews/PR284/man/formulations/index.html @@ -173,4 +173,4 @@ 0.0 0.02340899999999987 0.007743999999999858 - + diff --git a/previews/PR284/man/linearsolver/index.html b/previews/PR284/man/linearsolver/index.html index b2d89fa4..171d057c 100644 --- a/previews/PR284/man/linearsolver/index.html +++ b/previews/PR284/man/linearsolver/index.html @@ -1,3 +1,3 @@ Linear Solvers · ExaPF.jl

Linear Solver

Overview

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),
  • UMFPACK through the default Julia \ operator (CPU),
  • generic BiCGSTAB implementation [Vorst1992] (CPU/GPU),
  • or any linear solver wrapped in LinearAlgebra.

Preconditioning

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 block-Jacobi preconditioner used in ExaPF has been added to KrylovPreconditioners.jl

The Jacobian is partitioned into a dense block diagonal structure using Metis.jl, where each block is inverted to build our preconditioner P.

Dense block Jacobi preconditioner \\label{fig:preconditioner}

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.

  • Vorst1992Vorst, H. A. van der. 1992. “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 (2): 631–44
+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.

diff --git a/previews/PR284/man/powersystem/index.html b/previews/PR284/man/powersystem/index.html index f9a5f9c8..8d346632 100644 --- a/previews/PR284/man/powersystem/index.html +++ b/previews/PR284/man/powersystem/index.html @@ -16,4 +16,4 @@ 2 julia> PS.get(ps, PS.NumberOfSlackBuses()) -1 +1 diff --git a/previews/PR284/quickstart/index.html b/previews/PR284/quickstart/index.html index c512a6d2..96973332 100644 --- a/previews/PR284/quickstart/index.html +++ b/previews/PR284/quickstart/index.html @@ -15,14 +15,14 @@ #it 2: 5.88242e-01 #it 3: 4.88493e-03 #it 4: 1.39924e-06 -#it 5: 1.52178e-11 +#it 5: 1.57368e-11 Power flow has converged: true * #iterations: 5 - * Time Jacobian (s) ........: 0.4013 - * Time linear solver (s) ...: 0.0175 - * update (s) ............: 0.0171 - * ldiv (s) ..............: 0.0004 - * Time total (s) ...........: 0.8081

Implicitly, ExaPF has just proceed to the following operations:

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.3835
+  * Time linear solver (s) ...: 0.0017
+     * update (s) ............: 0.0015
+     * ldiv (s) ..............: 0.0002
+  * Time total (s) ...........: 0.7724

Implicitly, ExaPF has just proceed to the following operations:

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})) \,, \\ @@ -60,14 +60,14 @@ #it 2: 5.88242e-01 #it 3: 4.88493e-03 #it 4: 1.39924e-06 -#it 5: 1.52178e-11 +#it 5: 1.57368e-11 Power flow has converged: true * #iterations: 5 - * Time Jacobian (s) ........: 0.0045 - * Time linear solver (s) ...: 0.0173 - * update (s) ............: 0.0169 - * ldiv (s) ..............: 0.0004 - * Time total (s) ...........: 0.0222

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.0043
+  * Time linear solver (s) ...: 0.0017
+     * update (s) ............: 0.0015
+     * ldiv (s) ..............: 0.0002
+  * Time total (s) ...........: 0.0064

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
@@ -83,11 +83,11 @@
 #it 5: 9.96622e-12
 Power flow has converged: true
   * #iterations: 5
-  * Time Jacobian (s) ........: 2.2953
-  * Time linear solver (s) ...: 0.2142
+  * Time Jacobian (s) ........: 2.3728
+  * Time linear solver (s) ...: 0.2180
      * update (s) ............: 0.0000
-     * ldiv (s) ..............: 0.2142
-  * Time total (s) ...........: 3.2764

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 = BlockJacobiPreconditioner(jac_gpu, npartitions, CUDABackend());ERROR: UndefVarError: `BlockJacobiPreconditioner` not defined

You can attach the preconditioner to an BICGSTAB algorithm simply as

julia> linear_solver = ExaPF.KrylovBICGSTAB(jac_gpu; P=precond);ERROR: UndefVarError: `precond` not defined

(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
+     * ldiv (s) ..............: 0.2180
+  * Time total (s) ...........: 3.3683

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 = BlockJacobiPreconditioner(jac_gpu, npartitions, CUDABackend());ERROR: UndefVarError: `BlockJacobiPreconditioner` not defined

You can attach the preconditioner to an BICGSTAB algorithm simply as

julia> linear_solver = ExaPF.Bicgstab(jac_gpu; P=precond);ERROR: UndefVarError: `precond` not defined

(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
@@ -96,7 +96,7 @@
 Power flow has converged: true
   * #iterations: 5
   * Time Jacobian (s) ........: 0.0025
-  * Time linear solver (s) ...: 0.1146
+  * Time linear solver (s) ...: 0.1171
      * update (s) ............: 0.0000
-     * ldiv (s) ..............: 0.1146
-  * Time total (s) ...........: 0.1186
+ * ldiv (s) ..............: 0.1171 + * Time total (s) ...........: 0.1211 diff --git a/previews/PR284/search_index.js b/previews/PR284/search_index.js index b78e5840..4d328c35 100644 --- a/previews/PR284/search_index.js +++ b/previews/PR284/search_index.js @@ -1,3 +1,3 @@ var documenterSearchIndex = {"docs": -[{"location":"lib/formulations/","page":"Polar formulation","title":"Polar formulation","text":"CurrentModule = ExaPF\nconst PS = ExaPF.PowerSystem\nDocTestSetup = quote\n using ExaPF\nend","category":"page"},{"location":"lib/formulations/#Polar-formulation","page":"Polar formulation","title":"Polar formulation","text":"","category":"section"},{"location":"lib/formulations/#Generic-templates","page":"Polar formulation","title":"Generic templates","text":"","category":"section"},{"location":"lib/formulations/","page":"Polar formulation","title":"Polar formulation","text":"AbstractVariable\nAbstractFormulation\nState\nControl\n","category":"page"},{"location":"lib/formulations/#ExaPF.AbstractVariable","page":"Polar formulation","title":"ExaPF.AbstractVariable","text":"AbstractVariable\n\nVariables corresponding to a particular formulation.\n\n\n\n\n\n","category":"type"},{"location":"lib/formulations/#ExaPF.AbstractFormulation","page":"Polar formulation","title":"ExaPF.AbstractFormulation","text":"AbstractFormulation\n\nInterface between the data and the mathemical formulation.\n\n\n\n\n\n","category":"type"},{"location":"lib/formulations/#ExaPF.State","page":"Polar formulation","title":"ExaPF.State","text":"State <: AbstractVariable\n\nAll variables x depending on the variables Control u through the non-linear equation g(x u) = 0.\n\n\n\n\n\n","category":"type"},{"location":"lib/formulations/#ExaPF.Control","page":"Polar formulation","title":"ExaPF.Control","text":"Control <: AbstractVariable\n\nIndependent variables u used in the reduced-space formulation.\n\n\n\n\n\n","category":"type"},{"location":"lib/formulations/#Structure-and-variables","page":"Polar formulation","title":"Structure and variables","text":"","category":"section"},{"location":"lib/formulations/","page":"Polar formulation","title":"Polar formulation","text":"PolarForm\nBlockPolarForm\nload_polar\nNetworkStack\ninit!\n","category":"page"},{"location":"lib/formulations/#ExaPF.PolarForm","page":"Polar formulation","title":"ExaPF.PolarForm","text":"PolarForm{T, IT, VT, MT} <: AbstractPolarFormulation\n\nImplement the polar formulation associated to the network's equations.\n\nWrap a PS.PowerNetwork network to load the data on the target device (CPU() and CUDABackend() are currently supported).\n\nExample\n\njulia> const PS = ExaPF.PowerSystem;\n\njulia> network_data = PS.load_case(\"case9.m\");\n\njulia> polar = PolarForm(network_data, ExaPF.CPU())\nPolar formulation (instantiated on device CPU(false))\nNetwork characteristics:\n #buses: 9 (#slack: 1 #PV: 2 #PQ: 6)\n #generators: 3\n #lines: 9\ngiving a mathematical formulation with:\n #controls: 5\n #states : 14\n\n\n\n\n\n\n","category":"type"},{"location":"lib/formulations/#ExaPF.BlockPolarForm","page":"Polar formulation","title":"ExaPF.BlockPolarForm","text":"BlockPolarForm{T, IT, VT, MT} <: AbstractFormulation\n\nBlock polar formulation: duplicates k different polar models to evaluate them in parallel.\n\n\n\n\n\n","category":"type"},{"location":"lib/formulations/#ExaPF.load_polar","page":"Polar formulation","title":"ExaPF.load_polar","text":"load_polar(case, device=CPU(); dir=PS.EXADATA)\n\nLoad 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).\n\nExamples\n\njulia> polar = ExaPF.load_polar(\"case9\")\nPolar formulation (instantiated on device CPU(false))\nNetwork characteristics:\n #buses: 9 (#slack: 1 #PV: 2 #PQ: 6)\n #generators: 3\n #lines: 9\ngiving a mathematical formulation with:\n #controls: 5\n #states : 14\n\n\n\n\n\n\n","category":"function"},{"location":"lib/formulations/#ExaPF.NetworkStack","page":"Polar formulation","title":"ExaPF.NetworkStack","text":"NetworkStack{VT,VD,MT} <: AbstractNetworkStack{VT}\nNetworkStack(polar::PolarForm)\nNetworkStack(nbus::Int, ngen::Int, nlines::Int, VT::Type)\n\nStore the variables associated to the polar formulation. The variables are stored in the field input, ordered as follows\n\n input = [vmag ; vang ; pgen]\n\nThe object stores also intermediate variables needed in the expression tree, such as the LKMR basis ψ.\n\nNotes\n\nThe NetworkStack can be instantiated on the host or on the target device.\n\nExamples\n\njulia> polar = ExaPF.load_polar(\"case9\");\n\njulia> stack = ExaPF.NetworkStack(polar)\n21-elements NetworkStack{Vector{Float64}}\n\njulia> stack.vmag\n9-element Vector{Float64}:\n 1.0\n 1.0\n 1.0\n 1.0\n 1.0\n 1.0\n 1.0\n 1.0\n 1.0\n\n\n\n\n\n\n","category":"type"},{"location":"lib/formulations/#ExaPF.init!","page":"Polar formulation","title":"ExaPF.init!","text":"init!(polar::PolarForm, stack::NetworkStack)\n\nSet stack.input with the initial values specified in the base PS.PowerNetwork object.\n\n\n\n\n\n","category":"function"},{"location":"lib/formulations/","page":"Polar formulation","title":"Polar formulation","text":"The state and the control are defined as mapping:","category":"page"},{"location":"lib/formulations/","page":"Polar formulation","title":"Polar formulation","text":"mapping\n","category":"page"},{"location":"lib/formulations/#ExaPF.mapping","page":"Polar formulation","title":"ExaPF.mapping","text":"mapping(polar::PolarForm, ::Control)\n\nReturn the mapping associated to the Control() in NetworkStack according to the polar formulation PolarForm.\n\nExamples\n\njulia> polar = ExaPF.load_polar(\"case9\");\n\njulia> mapu = ExaPF.mapping(polar, Control())\n5-element Vector{Int64}:\n 1\n 2\n 3\n 20\n 21\n\n\n\n\n\n\nmapping(polar::PolarForm, ::State)\n\nReturn the mapping associated to the State() in NetworkStack according to the polar formulation PolarForm.\n\nExamples\n\njulia> polar = ExaPF.load_polar(\"case9\");\n\njulia> mapu = ExaPF.mapping(polar, State())\n14-element Vector{Int64}:\n 11\n 12\n 13\n 14\n 15\n 16\n 17\n 18\n 4\n 5\n 6\n 7\n 8\n 9\n\n\n\n\n\n\n","category":"function"},{"location":"lib/formulations/#Powerflow-solver","page":"Polar formulation","title":"Powerflow solver","text":"","category":"section"},{"location":"lib/formulations/","page":"Polar formulation","title":"Polar formulation","text":"run_pf\nnlsolve!\nNewtonRaphson\n","category":"page"},{"location":"lib/formulations/#ExaPF.run_pf","page":"Polar formulation","title":"ExaPF.run_pf","text":"run_pf(\n polar::PolarForm, stack::NetworkStack;\n rtol=1e-8, max_iter=20, verbose=0,\n)\n\nSolve 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.\n\nThe algorithm stops when a tolerance rtol or a maximum number of iterations maxiter is reached.\n\nArguments\n\npolar::AbstractFormulation: formulation of the power flow equation\nstack::NetworkStack: initial values in the network\n\nExamples\n\njulia> polar = ExaPF.load_polar(\"case9\");\n\njulia> stack = ExaPF.NetworkStack(polar);\n\njulia> conv = run_pf(polar, stack);\n\njulia> conv.has_converged\ntrue\n\njulia> conv.n_iterations\n4\n\n\n\n\n\n","category":"function"},{"location":"lib/formulations/#ExaPF.nlsolve!","page":"Polar formulation","title":"ExaPF.nlsolve!","text":"nlsolve!(\n algo::NewtonRaphson,\n jac::Jacobian,\n stack::NetworkStack;\n linear_solver=DirectSolver(jac.J),\n nl_buffer=NLBuffer(size(jac, 2)),\n)\n\nSolve the nonlinear system of equations g(x) = 0 with a NewtonRaphson algorithm. At each iteration, we update the variable x as\n\n x_k+1 = x_k - (g_k)^-1 g(x_k)\n\n\ntill g(x_k) ε_tol\n\nIn the implementation,\n\nthe function g is specified in jac.func,\nthe initial variable x_0 in stack::NetworkStack (with mapping jac.map),\nthe Jacobian g is computed automatically in jac, with automatic differentiation.\n\nNote that stack is modified inplace during the iterations of algorithm.\n\nThe Jacobian jac should be instantied before calling this function. By default, the linear system (g_k)^-1 g(x_k) is solved using a LU factorization. You can specify a different linear solver by changing the optional argument linear_solver.\n\nArguments\n\nalgo::NewtonRaphon: Newton-Raphson object, storing the options of the algorithm\njac::Jacobian: Stores the function g and its Jacobian g. The Jacobian is updated with automatic differentiation.\nstack::NetworkStack: initial values\nlinear_solver::AbstractLinearSolver: linear solver used to compute the Newton step\nnl_buffer::NLBuffer: buffer storing the residual vector and the descent direction Δx. Can be reused to avoid unecessary allocations.\n\nExamples\n\njulia> polar = ExaPF.load_polar(\"case9\");\n\njulia> powerflow = ExaPF.PowerFlowBalance(polar) ∘ ExaPF.PolarBasis(polar);\n\njulia> jx = ExaPF.Jacobian(polar, powerflow, State());\n\njulia> stack = ExaPF.NetworkStack(polar);\n\njulia> conv = ExaPF.nlsolve!(NewtonRaphson(), jx, stack);\n\njulia> conv.has_converged\ntrue\n\njulia> conv.n_iterations\n4\n\n\n\n\n\n","category":"function"},{"location":"lib/formulations/#ExaPF.NewtonRaphson","page":"Polar formulation","title":"ExaPF.NewtonRaphson","text":"NewtonRaphson <: AbstractNonLinearSolver\n\nNewton-Raphson algorithm.\n\nAttributes\n\nmaxiter::Int (default 20): maximum number of iterations\ntol::Float64 (default 1e-8): tolerance of the algorithm\nverbose::Int (default 0): verbosity level\n\n\n\n\n\n","category":"type"},{"location":"lib/formulations/#Constraints","page":"Polar formulation","title":"Constraints","text":"","category":"section"},{"location":"lib/formulations/","page":"Polar formulation","title":"Polar formulation","text":"The different parts of the polar formulation are implemented in the following AbstractExpression:","category":"page"},{"location":"lib/formulations/","page":"Polar formulation","title":"Polar formulation","text":"PolarBasis\nPowerFlowBalance\nVoltageMagnitudeBounds\nPowerGenerationBounds\nLineFlows\n","category":"page"},{"location":"lib/formulations/#ExaPF.PolarBasis","page":"Polar formulation","title":"ExaPF.PolarBasis","text":"PolarBasis{VI, MT} <: AbstractExpression\nPolarBasis(polar::AbstractPolarFormulation)\n\nImplement the LKMR nonlinear basis. Takes as input the voltage magnitudes vmag and the voltage angles vang and returns\n\n beginaligned\n psi_ell^C(v theta) = v^f v^t cos(theta_f - theta_t) quad forall ell = 1 cdots n_ell \n psi_ell^S(v theta) = v^f v^t sin(theta_f - theta_t) quad forall ell = 1 cdots n_ell \n psi_k(v theta) = v_k^2 quad forall k = 1 cdots n_b\n endaligned\n\nDimension: 2 * n_lines + n_bus\n\nComplexity\n\n3 n_lines + n_bus mul, n_lines cos and n_lines sin\n\nExamples\n\njulia> polar = ExaPF.load_polar(\"case9\");\n\njulia> stack = ExaPF.NetworkStack(polar);\n\njulia> basis = ExaPF.PolarBasis(polar)\nPolarBasis (AbstractExpression)\n\njulia> basis(stack)\n27-element Vector{Float64}:\n 1.0\n 1.0\n 1.0\n 1.0\n 1.0\n 1.0\n 1.0\n 1.0\n 1.0\n 0.0\n ⋮\n 1.0\n 1.0\n 1.0\n 1.0\n 1.0\n 1.0\n 1.0\n 1.0\n 1.0\n\n\n\n\n\n\n","category":"type"},{"location":"lib/formulations/#ExaPF.PowerFlowBalance","page":"Polar formulation","title":"ExaPF.PowerFlowBalance","text":"PowerFlowBalance{VT, MT}\nPowerFlowBalance(polar)\n\nImplement 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:\n\nbeginaligned\n p_i = v_i sum_j^n v_j (g_ijcos(theta_i - theta_j) + b_ijsin(theta_i - theta_j)) \n i PV PQ \n q_i = v_i sum_j^n v_j (g_ijsin(theta_i - theta_j) - b_ijcos(theta_i - theta_j)) \n i PQ\nendaligned\n\nRequire composition with PolarBasis.\n\nDimension: n_pv + 2 * n_pq\n\nComplexity\n\n2 SpMV\n\nExamples\n\njulia> polar = ExaPF.load_polar(\"case9\");\n\njulia> stack = ExaPF.NetworkStack(polar);\n\njulia> powerflow = ExaPF.PowerFlowBalance(polar) ∘ ExaPF.PolarBasis(polar);\n\njulia> round.(powerflow(stack); digits=6)\n14-element Vector{Float64}:\n -1.63\n -0.85\n 0.0\n 0.9\n 0.0\n 1.0\n 0.0\n 1.25\n -0.167\n 0.042\n -0.2835\n 0.171\n -0.2275\n 0.259\n\njulia> run_pf(polar, stack); # solve powerflow equations\n\njulia> isapprox(powerflow(stack), zeros(14); atol=1e-8)\ntrue\n\n\n\n\n\n\n","category":"type"},{"location":"lib/formulations/#ExaPF.VoltageMagnitudeBounds","page":"Polar formulation","title":"ExaPF.VoltageMagnitudeBounds","text":"VoltageMagnitudeBounds\n\nImplement 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:\n\nv_pq^ v_pq v_pq^ \n\nDimension: n_pq\n\nComplexity\n\n1 copyto\n\nNote\n\nIn the reduced space, the constraints on the voltage magnitudes at PV nodes v_pv are taken into account when bounding the control u.\n\nExamples\n\njulia> polar = ExaPF.load_polar(\"case9\");\n\njulia> stack = ExaPF.NetworkStack(polar);\n\njulia> voltage_pq = ExaPF.VoltageMagnitudeBounds(polar);\n\njulia> voltage_pq(stack)\n6-element Vector{Float64}:\n 1.0\n 1.0\n 1.0\n 1.0\n 1.0\n 1.0\n\n\n\n\n\n\n","category":"type"},{"location":"lib/formulations/#ExaPF.PowerGenerationBounds","page":"Polar formulation","title":"ExaPF.PowerGenerationBounds","text":"PowerGenerationBounds{VT, MT}\nPowerGenerationBounds(polar)\n\nConstraints 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\n\np_gref^ p_gref p_gref^ \nC_g q_g^ C_g q_g C_g q_g^ \n\nRequire composition with PolarBasis.\n\nDimension: n_pv + 2 n_ref\n\nComplexity\n\n1 copyto, 1 SpMV\n\nExamples\n\njulia> polar = ExaPF.load_polar(\"case9\");\n\njulia> stack = ExaPF.NetworkStack(polar);\n\njulia> run_pf(polar, stack); # solve powerflow equations\n\njulia> power_generators = ExaPF.PowerGenerationBounds(polar) ∘ ExaPF.PolarBasis(polar);\n\njulia> round.(power_generators(stack); digits=6)\n4-element Vector{Float64}:\n 0.719547\n 0.24069\n 0.144601\n -0.03649\n\n\n\n\n\n\n","category":"type"},{"location":"lib/formulations/#ExaPF.LineFlows","page":"Polar formulation","title":"ExaPF.LineFlows","text":"LineFlows{VT, MT}\nLineFlows(polar)\n\nImplement thermal limit constraints on the lines of the network.\n\nRequire composition with PolarBasis.\n\nDimension: 2 * n_lines\n\nComplexity\n\n4 SpMV, 4 * n_lines quadratic, 2 * n_lines add\n\nExamples\n\njulia> polar = ExaPF.load_polar(\"case9\");\n\njulia> stack = ExaPF.NetworkStack(polar);\n\njulia> run_pf(polar, stack); # solve powerflow equations\n\njulia> line_flows = ExaPF.LineFlows(polar) ∘ ExaPF.PolarBasis(polar);\n\njulia> round.(line_flows(stack); digits=6)\n18-element Vector{Float64}:\n 0.575679\n 0.094457\n 0.379983\n 0.723832\n 0.060169\n 0.588673\n 2.657418\n 0.748943\n 0.295351\n 0.560817\n 0.112095\n 0.38625\n 0.728726\n 0.117191\n 0.585164\n 2.67781\n 0.726668\n 0.215497\n\n\n\n\n\n\n","category":"type"},{"location":"lib/formulations/#Objective","page":"Polar formulation","title":"Objective","text":"","category":"section"},{"location":"lib/formulations/","page":"Polar formulation","title":"Polar formulation","text":"The production costs is given in the AbstractExpression CostFunction:","category":"page"},{"location":"lib/formulations/","page":"Polar formulation","title":"Polar formulation","text":"CostFunction","category":"page"},{"location":"lib/formulations/#ExaPF.CostFunction","page":"Polar formulation","title":"ExaPF.CostFunction","text":"CostFunction{VT, MT} <: AutoDiff.AbstractExpression\nCostFunction(polar)\n\nImplement the quadratic cost function for OPF\n\n _g=1^n_g c_2g p_g^2 + c_1g p_g + c_0g\n\nRequire composition with PolarBasis to evaluate the cost of the reference generator.\n\nDimension: 1\n\nComplexity\n\n1 SpMV, 1 sum\n\nExamples\n\njulia> polar = ExaPF.load_polar(\"case9\");\n\njulia> stack = ExaPF.NetworkStack(polar);\n\njulia> cost = ExaPF.CostFunction(polar) ∘ ExaPF.PolarBasis(polar);\n\njulia> cost(stack)\n1-element Vector{Float64}:\n 4509.0275\n\n\n\n\n\n\n","category":"type"},{"location":"lib/formulations/#Composition-of-expressions","page":"Polar formulation","title":"Composition of expressions","text":"","category":"section"},{"location":"lib/formulations/","page":"Polar formulation","title":"Polar formulation","text":"The different expressions can be combined together in several different ways.","category":"page"},{"location":"lib/formulations/","page":"Polar formulation","title":"Polar formulation","text":"MultiExpressions\nComposedExpressions","category":"page"},{"location":"lib/formulations/#ExaPF.MultiExpressions","page":"Polar formulation","title":"ExaPF.MultiExpressions","text":"MultiExpressions <: AbstractExpression\n\nImplement expressions concatenation. Takes as input a vector of expressions [expr1,...,exprN] and concatenate them in a single expression mexpr, such that\n\n mexpr(x) = [expr1(x) ; expr2(x) ; ... ; exprN(x)]\n\n\n\n\n\n\n","category":"type"},{"location":"lib/formulations/#ExaPF.ComposedExpressions","page":"Polar formulation","title":"ExaPF.ComposedExpressions","text":"ComposedExpressions{Expr1<:PolarBasis, Expr2} <: AbstractExpression\n\nImplement expression composition. Takes as input two expressions expr1 and expr2 and returns a composed expression cexpr such that ``` cexpr(x) = expr2 ∘ expr1(x)\n\nNotes\n\nCurrently, only PolarBasis is supported for expr1.\n\n\n\n\n\n","category":"type"},{"location":"lib/autodiff/","page":"AutoDiff","title":"AutoDiff","text":"CurrentModule = ExaPF.AutoDiff","category":"page"},{"location":"lib/autodiff/#AutoDiffAPI","page":"AutoDiff","title":"AutoDiff","text":"","category":"section"},{"location":"lib/autodiff/#Variables","page":"AutoDiff","title":"Variables","text":"","category":"section"},{"location":"lib/autodiff/","page":"AutoDiff","title":"AutoDiff","text":"AbstractStack\n","category":"page"},{"location":"lib/autodiff/#ExaPF.AutoDiff.AbstractStack","page":"AutoDiff","title":"ExaPF.AutoDiff.AbstractStack","text":"AbstractStack{VT}\n\nAbstract variable storing the inputs and the intermediate values in the expression tree.\n\n\n\n\n\n","category":"type"},{"location":"lib/autodiff/#Expressions","page":"AutoDiff","title":"Expressions","text":"","category":"section"},{"location":"lib/autodiff/","page":"AutoDiff","title":"AutoDiff","text":"AbstractExpression\nadjoint!\n","category":"page"},{"location":"lib/autodiff/#ExaPF.AutoDiff.AbstractExpression","page":"AutoDiff","title":"ExaPF.AutoDiff.AbstractExpression","text":"AbstractExpression\n\nAbstract 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).\n\nForward mode\n\nThe direct evaluation of the function f is implemented as\n\n(expr::AbstractExpression)(output::VT, stack::AbstractStack{VT}) where VT<:AbstractArray\n\n\nthe input being specified in stack, the results being stored in the array output.\n\nReverse mode\n\nThe adjoint of the function is specified by the function adjoint!, with the signature:\n\nadjoint!(expr::AbstractExpression, ∂stack::AbstractStack{VT}, stack::AbstractStack{VT}, ̄v::VT) where VT<:AbstractArray\n\n\nThe variable stack stores the result of the direct evaluation, and is not modified in adjoint!. The results are stored inside the adjoint stack ∂stack.\n\n\n\n\n\n","category":"type"},{"location":"lib/autodiff/#ExaPF.AutoDiff.adjoint!","page":"AutoDiff","title":"ExaPF.AutoDiff.adjoint!","text":"adjoint!(expr::AbstractExpression, ∂stack::AbstractStack{VT}, stack::AbstractStack{VT}, ̄v::VT) where VT<:AbstractArray\n\nCompute 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!.\n\n\n\n\n\n","category":"function"},{"location":"lib/autodiff/#First-and-second-order-derivatives","page":"AutoDiff","title":"First and second-order derivatives","text":"","category":"section"},{"location":"lib/autodiff/","page":"AutoDiff","title":"AutoDiff","text":"AbstractJacobian\nAbstractHessianProd\nAbstractFullHessian","category":"page"},{"location":"lib/autodiff/#ExaPF.AutoDiff.AbstractJacobian","page":"AutoDiff","title":"ExaPF.AutoDiff.AbstractJacobian","text":"AbstractJacobian\n\nAutomatic differentiation for the compressed Jacobian of any nonlinear constraint h(x).\n\n\n\n\n\n","category":"type"},{"location":"lib/autodiff/#ExaPF.AutoDiff.AbstractHessianProd","page":"AutoDiff","title":"ExaPF.AutoDiff.AbstractHessianProd","text":"AbstractHessianProd\n\nReturns the adjoint-Hessian-vector product λ^ H v of any nonlinear constraint h(x).\n\n\n\n\n\n","category":"type"},{"location":"lib/autodiff/#ExaPF.AutoDiff.AbstractFullHessian","page":"AutoDiff","title":"ExaPF.AutoDiff.AbstractFullHessian","text":"AbstractHessianProd\n\nFull sparse Hessian H of any nonlinear constraint h(x).\n\n\n\n\n\n","category":"type"},{"location":"lib/autodiff/#Utils","page":"AutoDiff","title":"Utils","text":"","category":"section"},{"location":"lib/autodiff/","page":"AutoDiff","title":"AutoDiff","text":"seed!\nseed_coloring!\npartials!\nset_value!","category":"page"},{"location":"lib/autodiff/#ExaPF.AutoDiff.seed!","page":"AutoDiff","title":"ExaPF.AutoDiff.seed!","text":"seed!(\n H::AbstractHessianProd,\n v::AbstractVector{T},\n) where {T}\n\nSeed the duals with v to compute the Hessian vector product λ^ H v.\n\n\n\n\n\n","category":"function"},{"location":"lib/autodiff/#ExaPF.AutoDiff.seed_coloring!","page":"AutoDiff","title":"ExaPF.AutoDiff.seed_coloring!","text":"seed_coloring!(\n M::Union{AbstractJacobian, AbstractFullHessian}\n coloring::AbstractVector,\n)\n\nSeed the duals with the coloring based seeds to compute the Jacobian or Hessian M.\n\n\n\n\n\n","category":"function"},{"location":"lib/autodiff/#ExaPF.AutoDiff.partials!","page":"AutoDiff","title":"ExaPF.AutoDiff.partials!","text":"partials!(jac::AbstractJacobian)\n\nExtract partials from Jacobian jac in jac.J.\n\n\n\n\n\npartials!(hess::AbstractFullHessian)\n\nExtract partials from Hessian hess into hess.H.\n\n\n\n\n\n","category":"function"},{"location":"lib/autodiff/#ExaPF.AutoDiff.set_value!","page":"AutoDiff","title":"ExaPF.AutoDiff.set_value!","text":"set_value!(\n jac,\n primals::AbstractVector{T}\n) where {T}\n\nSet values of ForwardDiff.Dual numbers in jac to primals.\n\n\n\n\n\n","category":"function"},{"location":"man/powersystem/","page":"PowerSystem","title":"PowerSystem","text":"CurrentModule = ExaPF\nDocTestSetup = quote\n using ExaPF\n const PS = ExaPF.PowerSystem\nend\nDocTestFilters = [r\"ExaPF\"]","category":"page"},{"location":"man/powersystem/#PowerSystem","page":"PowerSystem","title":"PowerSystem","text":"","category":"section"},{"location":"man/powersystem/","page":"PowerSystem","title":"PowerSystem","text":"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.","category":"page"},{"location":"man/powersystem/","page":"PowerSystem","title":"PowerSystem","text":"We usually load the PowerSystem system submodule with the alias PS:","category":"page"},{"location":"man/powersystem/","page":"PowerSystem","title":"PowerSystem","text":"julia> PS = ExaPF.PowerSystem\n","category":"page"},{"location":"man/powersystem/#Description","page":"PowerSystem","title":"Description","text":"","category":"section"},{"location":"man/powersystem/","page":"PowerSystem","title":"PowerSystem","text":"The electrical power system is represented as a linear, lumped network which has to satisfy the Kirchhoff laws:","category":"page"},{"location":"man/powersystem/","page":"PowerSystem","title":"PowerSystem","text":" bmi = bmYbmv ","category":"page"},{"location":"man/powersystem/","page":"PowerSystem","title":"PowerSystem","text":"where bmi bmv in mathbbC^N_B are the current and voltage vectors associated to the system and bmY in mathbbC^N_B times N_B is the admittance matrix. These equations are often rewritten in terms of apparent powers:","category":"page"},{"location":"man/powersystem/","page":"PowerSystem","title":"PowerSystem","text":" bms = bmp + jbmq = textitdiag(bmv^*) bmYbmv","category":"page"},{"location":"man/powersystem/","page":"PowerSystem","title":"PowerSystem","text":"Using polar representation of the voltage vector, such as bmv = ve^j theta, each bus i=1 cdots N_B must satisfy the power balance equations:","category":"page"},{"location":"man/powersystem/","page":"PowerSystem","title":"PowerSystem","text":"beginaligned\n p_i = v_i sum_j^n v_j (g_ijcos(theta_i - theta_j) + b_ijsin(theta_i - theta_j)) \n q_i = v_i sum_j^n v_j (g_ijsin(theta_i - theta_j) - b_ijcos(theta_i - theta_j)) \nendaligned","category":"page"},{"location":"man/powersystem/","page":"PowerSystem","title":"PowerSystem","text":"where each bus i has variables p_i q_i v_i theta_i and the topology of the network is defined by a non-negative value of the admittance between two buses i and j, y_ij = g_ij + ib_ij.","category":"page"},{"location":"man/powersystem/#The-PowerNetwork-Object","page":"PowerSystem","title":"The PowerNetwork Object","text":"","category":"section"},{"location":"man/powersystem/","page":"PowerSystem","title":"PowerSystem","text":"Currently we can create a PS.PowerNetwork object by parsing a MATPOWER data file.","category":"page"},{"location":"man/powersystem/","page":"PowerSystem","title":"PowerSystem","text":"julia> datafile = \"case9.m\";\n\njulia> ps = PS.load_case(datafile)\nPowerNetwork object with:\n Buses: 9 (Slack: 1. PV: 2. PQ: 6)\n Generators: 3.\n","category":"page"},{"location":"man/powersystem/","page":"PowerSystem","title":"PowerSystem","text":"Then, using multiple dispatch, we have defined a set of abstract data types and getter functions which allow us to retrieve information from the PowerNetwork object","category":"page"},{"location":"man/powersystem/","page":"PowerSystem","title":"PowerSystem","text":"julia> PS.get(ps, PS.NumberOfPQBuses())\n6\n\njulia> PS.get(ps, PS.NumberOfPVBuses())\n2\n\njulia> PS.get(ps, PS.NumberOfSlackBuses())\n1","category":"page"},{"location":"man/linearsolver/","page":"Linear Solvers","title":"Linear Solvers","text":"CurrentModule = ExaPF\nDocTestSetup = quote\n using ExaPF\n const Precondition = ExaPF.Precondition\n const Iterative = ExaPF.Iterative\nend\nDocTestFilters = [r\"ExaPF\"]","category":"page"},{"location":"man/linearsolver/#Linear-Solver","page":"Linear Solvers","title":"Linear Solver","text":"","category":"section"},{"location":"man/linearsolver/#Overview","page":"Linear Solvers","title":"Overview","text":"","category":"section"},{"location":"man/linearsolver/","page":"Linear Solvers","title":"Linear Solvers","text":"As mentioned before, a linear solver is required to compute the Newton step in","category":"page"},{"location":"man/linearsolver/","page":"Linear Solvers","title":"Linear Solvers","text":"dx .= jacobian(x)\\f(x)","category":"page"},{"location":"man/linearsolver/","page":"Linear Solvers","title":"Linear Solvers","text":"Our package supports the following linear solvers:","category":"page"},{"location":"man/linearsolver/","page":"Linear Solvers","title":"Linear Solvers","text":"cuSOLVER with csrlsvqr (GPU),\nKrylov.jl with dqgmres and bicgstab (CPU/GPU),\nUMFPACK through the default Julia \\ operator (CPU),\ngeneric BiCGSTAB implementation [Vorst1992] (CPU/GPU),\nor any linear solver wrapped in LinearAlgebra.","category":"page"},{"location":"man/linearsolver/#Preconditioning","page":"Linear Solvers","title":"Preconditioning","text":"","category":"section"},{"location":"man/linearsolver/","page":"Linear Solvers","title":"Linear Solvers","text":"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.","category":"page"},{"location":"man/linearsolver/","page":"Linear Solvers","title":"Linear Solvers","text":"The block-Jacobi preconditioner used in ExaPF has been added to KrylovPreconditioners.jl","category":"page"},{"location":"man/linearsolver/","page":"Linear Solvers","title":"Linear Solvers","text":"The Jacobian is partitioned into a dense block diagonal structure using Metis.jl, where each block is inverted to build our preconditioner P.","category":"page"},{"location":"man/linearsolver/","page":"Linear Solvers","title":"Linear Solvers","text":"(Image: Dense block Jacobi preconditioner \\label{fig:preconditioner})","category":"page"},{"location":"man/linearsolver/","page":"Linear Solvers","title":"Linear Solvers","text":"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.","category":"page"},{"location":"man/linearsolver/","page":"Linear Solvers","title":"Linear Solvers","text":"CUDA.@sync pivot, info = CUDA.CUBLAS.getrf_batched!(blocks, true)\nCUDA.@sync pivot, info, p.cuJs = CUDA.CUBLAS.getri_batched(blocks, pivot)","category":"page"},{"location":"man/linearsolver/","page":"Linear Solvers","title":"Linear Solvers","text":"Assuming that other vendors will provide such batched BLAS APIs, this code is portable to other GPU architectures.","category":"page"},{"location":"man/linearsolver/","page":"Linear Solvers","title":"Linear Solvers","text":"[Vorst1992]: Vorst, H. A. van der. 1992. “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 (2): 631–44","category":"page"},{"location":"artifact/#ExaData-Artifact","page":"ExaData Artifact","title":"ExaData Artifact","text":"","category":"section"},{"location":"artifact/","page":"ExaData Artifact","title":"ExaData Artifact","text":"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.","category":"page"},{"location":"tutorials/direct_solver/","page":"Power flow: direct solver","title":"Power flow: direct solver","text":"CurrentModule = ExaPF\nDocTestSetup = quote\n using ExaPF\n const LS = ExaPF.LinearSolvers\nend\nDocTestFilters = [r\"ExaPF\"]","category":"page"},{"location":"tutorials/direct_solver/","page":"Power flow: direct solver","title":"Power flow: direct solver","text":"using ExaPF\nusing KLU\nusing LinearAlgebra\nconst LS = ExaPF.LinearSolvers\n","category":"page"},{"location":"tutorials/direct_solver/#Direct-solvers-for-power-flow","page":"Power flow: direct solver","title":"Direct solvers for power flow","text":"","category":"section"},{"location":"tutorials/direct_solver/","page":"Power flow: direct solver","title":"Power flow: direct solver","text":"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","category":"page"},{"location":"tutorials/direct_solver/","page":"Power flow: direct solver","title":"Power flow: direct solver","text":"g(x p) = 0","category":"page"},{"location":"tutorials/direct_solver/","page":"Power flow: direct solver","title":"Power flow: direct solver","text":"where g mathbbR^n_x times mathbbR^n_p to mathbbR^n_x is a nonlinear function encoding the power flow equations.","category":"page"},{"location":"tutorials/direct_solver/","page":"Power flow: direct solver","title":"Power flow: direct solver","text":"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","category":"page"},{"location":"tutorials/direct_solver/","page":"Power flow: direct solver","title":"Power flow: direct solver","text":"(nabla_x g_k) Delta x = - g_k","category":"page"},{"location":"tutorials/direct_solver/","page":"Power flow: direct solver","title":"Power flow: direct solver","text":"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.","category":"page"},{"location":"tutorials/direct_solver/","page":"Power flow: direct solver","title":"Power flow: direct solver","text":"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","category":"page"},{"location":"tutorials/direct_solver/","page":"Power flow: direct solver","title":"Power flow: direct solver","text":"P (nabla_x g) Q = LU","category":"page"},{"location":"tutorials/direct_solver/","page":"Power flow: direct solver","title":"Power flow: direct solver","text":"Once the matrix factorized, solving the linear system just translates to perform two backsolves with the triangular matrices L and U.","category":"page"},{"location":"tutorials/direct_solver/","page":"Power flow: direct solver","title":"Power flow: direct solver","text":"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.","category":"page"},{"location":"tutorials/direct_solver/#UMFPACK-(CPU,-default)","page":"Power flow: direct solver","title":"UMFPACK (CPU, default)","text":"","category":"section"},{"location":"tutorials/direct_solver/","page":"Power flow: direct solver","title":"Power flow: direct solver","text":"By default, ExaPF employs the linear solver UMFPACK to solve the linear system, as UMFPACK is shipped automatically in Julia.","category":"page"},{"location":"tutorials/direct_solver/","page":"Power flow: direct solver","title":"Power flow: direct solver","text":"In the LinearSolvers submodule, this is how the wrapper is implemented:","category":"page"},{"location":"tutorials/direct_solver/","page":"Power flow: direct solver","title":"Power flow: direct solver","text":"struct DirectSolver{Fac} <: AbstractLinearSolver\n factorization::Fac\nend\nDirectSolver(J::AbstractMatrix) = DirectSolver(lu(J))\n","category":"page"},{"location":"tutorials/direct_solver/","page":"Power flow: direct solver","title":"Power flow: direct solver","text":"By default, the constructor takes as input the initial Jacobian J and factorizes it by calling lu(J), which in Julia translates to a factorization with UMFPACK. Then, inside the function nlsolve! we refactorize the matrix at each iteration by calling the function LinearSolvers.update!","category":"page"},{"location":"tutorials/direct_solver/","page":"Power flow: direct solver","title":"Power flow: direct solver","text":"function update!(s::DirectSolver, J::AbstractMatrix)\n LinearAlgebra.lu!(s.factorization, J) # Update factorization inplace\nend","category":"page"},{"location":"tutorials/direct_solver/","page":"Power flow: direct solver","title":"Power flow: direct solver","text":"This function uses the function LinearAlgebra.lu! to update the factorization inplace. The backsolve is computed by calling the LinearAlgebra.ldiv! function:","category":"page"},{"location":"tutorials/direct_solver/","page":"Power flow: direct solver","title":"Power flow: direct solver","text":"function ldiv!(s::DirectSolver, y::AbstractVector, J::AbstractMatrix, x::AbstractVector)\n LinearAlgebra.ldiv!(y, s.factorization, x) # Forward-backward solve\n return 0\nend","category":"page"},{"location":"tutorials/direct_solver/","page":"Power flow: direct solver","title":"Power flow: direct solver","text":"We notice that the code has been designed to support any factorization routines implementing the two routines LinearAlgebra.lu! and LinearAlgebra.ldiv!.","category":"page"},{"location":"tutorials/direct_solver/","page":"Power flow: direct solver","title":"Power flow: direct solver","text":"Before comparing with other linear solvers, we solve a large scale power flow instance with UMFPACK to give us a reference.","category":"page"},{"location":"tutorials/direct_solver/","page":"Power flow: direct solver","title":"Power flow: direct solver","text":"polar = ExaPF.load_polar(\"case9241pegase.m\")\nstack = ExaPF.NetworkStack(polar)\npf_solver = NewtonRaphson(tol=1e-10, verbose=2) # power flow solver\nfunc = ExaPF.PowerFlowBalance(polar) ∘ ExaPF.PolarBasis(polar) # power flow func\njx = ExaPF.Jacobian(polar, func, State()) # init AD\nExaPF.nlsolve!(pf_solver, jx, stack)","category":"page"},{"location":"tutorials/direct_solver/#KLU-(CPU)","page":"Power flow: direct solver","title":"KLU (CPU)","text":"","category":"section"},{"location":"tutorials/direct_solver/","page":"Power flow: direct solver","title":"Power flow: direct solver","text":"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!","category":"page"},{"location":"tutorials/direct_solver/","page":"Power flow: direct solver","title":"Power flow: direct solver","text":"LinearAlgebra.lu!(K::KLU.KLUFactorization, J) = KLU.klu!(K, J)\n","category":"page"},{"location":"tutorials/direct_solver/","page":"Power flow: direct solver","title":"Power flow: direct solver","text":"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:","category":"page"},{"location":"tutorials/direct_solver/","page":"Power flow: direct solver","title":"Power flow: direct solver","text":"klu_factorization = KLU.klu(jx.J)\nklu_solver = LS.DirectSolver(klu_factorization)\n","category":"page"},{"location":"tutorials/direct_solver/","page":"Power flow: direct solver","title":"Power flow: direct solver","text":"and pass it to nlsolve!:","category":"page"},{"location":"tutorials/direct_solver/","page":"Power flow: direct solver","title":"Power flow: direct solver","text":"ExaPF.init!(polar, stack) # reinit stack\nExaPF.nlsolve!(pf_solver, jx, stack; linear_solver=klu_solver)\n","category":"page"},{"location":"tutorials/direct_solver/","page":"Power flow: direct solver","title":"Power flow: direct solver","text":"We observe KLU reduces considerably the time spent in the linear solver.","category":"page"},{"location":"tutorials/direct_solver/#cusolverRF-(CUDA)","page":"Power flow: direct solver","title":"cusolverRF (CUDA)","text":"","category":"section"},{"location":"tutorials/direct_solver/","page":"Power flow: direct solver","title":"Power flow: direct solver","text":"cusolverRF is an efficient LU refactorization routine implemented in CUDA. It is wrapped in Julia inside the package CUSOLVERRF.jl:","category":"page"},{"location":"tutorials/direct_solver/","page":"Power flow: direct solver","title":"Power flow: direct solver","text":"using CUSOLVERRF","category":"page"},{"location":"tutorials/direct_solver/","page":"Power flow: direct solver","title":"Power flow: direct solver","text":"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.","category":"page"},{"location":"tutorials/direct_solver/","page":"Power flow: direct solver","title":"Power flow: direct solver","text":"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.","category":"page"},{"location":"tutorials/direct_solver/","page":"Power flow: direct solver","title":"Power flow: direct solver","text":"CUSOLVERRF.jl follows the LinearAlgebra's interface, so we can use it directly in ExaPF. We first have to instantiate everything on the GPU:","category":"page"},{"location":"tutorials/direct_solver/","page":"Power flow: direct solver","title":"Power flow: direct solver","text":"using CUDA\npolar_gpu = ExaPF.load_polar(\"case9241pegase.m\", CUDABackend())\nstack_gpu = ExaPF.NetworkStack(polar_gpu)\nfunc_gpu = ExaPF.PowerFlowBalance(polar_gpu) ∘ ExaPF.PolarBasis(polar_gpu)\njx_gpu = ExaPF.Jacobian(polar_gpu, func_gpu, State()) # init AD","category":"page"},{"location":"tutorials/direct_solver/","page":"Power flow: direct solver","title":"Power flow: direct solver","text":"We can instantiate a new cusolverRF's instance as","category":"page"},{"location":"tutorials/direct_solver/","page":"Power flow: direct solver","title":"Power flow: direct solver","text":"rf_fac = CUSOLVERRF.RFLU(jx_gpu.J)\nrf_solver = LS.DirectSolver(rf_fac)\n","category":"page"},{"location":"tutorials/direct_solver/","page":"Power flow: direct solver","title":"Power flow: direct solver","text":"Then, we are able to solve the power flow entirely on the GPU, simply as","category":"page"},{"location":"tutorials/direct_solver/","page":"Power flow: direct solver","title":"Power flow: direct solver","text":"ExaPF.nlsolve!(pf_solver, jx_gpu, stack_gpu; linear_solver=rf_solver)\n","category":"page"},{"location":"man/benchmark/","page":"Benchmark","title":"Benchmark","text":"CurrentModule = ExaPF\nDocTestSetup = quote\n using ExaPF\nend\nDocTestFilters = [r\"ExaPF\"]","category":"page"},{"location":"man/benchmark/#Benchmark","page":"Benchmark","title":"Benchmark","text":"","category":"section"},{"location":"man/benchmark/","page":"Benchmark","title":"Benchmark","text":"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.","category":"page"},{"location":"man/benchmark/","page":"Benchmark","title":"Benchmark","text":"Running","category":"page"},{"location":"man/benchmark/","page":"Benchmark","title":"Benchmark","text":"julia --project benchmark/benchmarks.jl KrylovBICGSTAB CUDABackend case300.m","category":"page"},{"location":"man/benchmark/","page":"Benchmark","title":"Benchmark","text":"yields","category":"page"},{"location":"man/benchmark/","page":"Benchmark","title":"Benchmark","text":"KrylovBICGSTAB, CUDABackend, case300.m, 69.0, 3.57, 43.7, true","category":"page"},{"location":"man/benchmark/","page":"Benchmark","title":"Benchmark","text":"The first three fields are the settings of the benchmark run. They are followed by three timings in milliseconds:","category":"page"},{"location":"man/benchmark/","page":"Benchmark","title":"Benchmark","text":"the time taken by the Newton-Raphson algorithm to solve the power flow,\nthe timings for the Jacobian accumulation using AutoDiff,\nand the time for the linear solver, including the preconditioner.","category":"page"},{"location":"man/benchmark/","page":"Benchmark","title":"Benchmark","text":"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","category":"page"},{"location":"man/benchmark/","page":"Benchmark","title":"Benchmark","text":"powerflow_solver = NewtonRaphson(tol=ntol)","category":"page"},{"location":"man/benchmark/","page":"Benchmark","title":"Benchmark","text":"to one of the following options:","category":"page"},{"location":"man/benchmark/","page":"Benchmark","title":"Benchmark","text":"powerflow_solver = NewtonRaphson(tol=ntol, verbose=VERBOSE_LEVEL_NONE)\npowerflow_solver = NewtonRaphson(tol=ntol, verbose=VERBOSE_LEVEL_LOW)\npowerflow_solver = NewtonRaphson(tol=ntol, verbose=VERBOSE_LEVEL_MEDIUM)\npowerflow_solver = NewtonRaphson(tol=ntol, verbose=VERBOSE_LEVEL_HIGH)","category":"page"},{"location":"man/benchmark/","page":"Benchmark","title":"Benchmark","text":"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.","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"CurrentModule = ExaPF\nDocTestSetup = quote\n using ExaPF\n using LazyArtifacts\n import ExaPF: AutoDiff\nend","category":"page"},{"location":"man/formulations/#Formulations","page":"Polar formulation","title":"Formulations","text":"","category":"section"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"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.","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"In what follows, we denote by v in mathbbR^n_b the voltage magnitudes, theta in mathbbR^n_b the voltage angles and p_g q_g in mathbbR^n_g the active and reactive power generations. The active and reactive loads are denoted respectively by p_d q_d in mathbbR^n_b.","category":"page"},{"location":"man/formulations/#Power-flow-model","page":"Polar formulation","title":"Power flow model","text":"","category":"section"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"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 mathbbR^n_b times mathbbR^n_b to mathbbR^2n_ell + n_b. If we introduce the intermediate expressions","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":" psi_ell^C(v theta) = v^f v^t cos(theta_f - theta_t) quad forall ell = 1 cdots n_ell \n psi_ell^S(v theta) = v^f v^t sin(theta_f - theta_t) quad forall ell = 1 cdots n_ell \n psi_k(v theta) = v_k^2 quad forall k = 1 cdots n_b","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"the basis psi is defined as","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":" psi(v theta) = psi_ell^C(v theta)^top psi_ell^S(v theta)^top psi_k(v theta)^top ","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"The basis psi encodes all the nonlinearities in the problem. For instance, the power flow equations rewrite directly as","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":" beginbmatrix\n C_g p_g - p_d \n C_g q_g - q_d\n endbmatrix\n +\n underbrace\n beginbmatrix\n - hatG^c - hatB^s -G^d \n hatB^c - hatG^s B^d\n endbmatrix\n _M\n psi(v theta)\n = 0","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"with C_g in mathbbR^n_b times n_g the bus-generators incidence matrix, and the matrices B G defined from the admittance matrix Y_b of the network.","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"Similarly, the line flows rewrite","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":" beginbmatrix\n s_p^f s_q^f\n endbmatrix\n =\n overbrace\n beginbmatrix\n G_ft B_ft G_ff C_f^top \n -B_ft G_ft -B_ff C_f^top\n endbmatrix\n ^L_line^f\n psi(v theta) \n beginbmatrix\n s_p^t s_q^t\n endbmatrix\n =\n underbrace\n beginbmatrix\n G_tf B_tf G_tt C_t^top \n -B_tf G_tf -B_tt C_t^top\n endbmatrix\n _L_line^t\n psi(v theta)","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"with C_f in mathbbR^n_b times n_ell the bus-from incidence matrix and C_t in mathbbR^n_b times n_ell the bus-to incidence matrix. Then, the line flows constraints write directly with the quadratic expressions:","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":" (s_p^f)^2 + (s_q^f)^2 leq (s^max)^2 quad \n (s_p^t)^2 + (s_q^t)^2 leq (s^max)^2 quad ","category":"page"},{"location":"man/formulations/#Why-is-this-model-advantageous?","page":"Polar formulation","title":"Why is this model advantageous?","text":"","category":"section"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"Implementing the model LTMR2020 is not difficult once the basis function psi has been defined. Indeed, if we select a subset of the power flow equations (as usual, associated to the active injections at PV nodes, and active and reactive injections at PQ nodes), we get","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":" C_eq p_g + M_eq psi + tau = 0","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"with C_eq defined from the bus-generator incidence matrix C_g, M_eq a subset of the matrix M, tau a constant depending on the loads in the problem. Note that C_eq and M_eq are sparse matrices, so the expression can be implemented efficiently with sparse linear algebra operations (2 SpMV operations, 2 vector additions). The same holds true for the line flow constraints, evaluated with 2 SpMV operations:","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":" s^f = L_line^f psi quad\n s^t = L_line^t psi ","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"In ExaPF, all nonlinear expressions are written as linear operations depending on the nonlinear basis psi. By doing so, all the unstructured sparsity of the power flow problem is directly handled inside the sparse linear algebra library (cusparse on CUDA GPU, SuiteSparse on the CPU).","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"In what follows, we detail the implementation of the LTMR2020 model in ExaPF.","category":"page"},{"location":"man/formulations/#How-to-instantiate-the-inputs?","page":"Polar formulation","title":"How to instantiate the inputs?","text":"","category":"section"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"We have implemented the LTMR2020 model in ExaPF, both on the CPU and on CUDA GPU. All the operations have been rewritten in a vectorized fashion. Every model depends on inputs we propagate forward with functions. In ExaPF, the inputs will be specified in a NetworkStack <: AbstractStack. The functions will be implemented as AutoDiff.AbstractExpression.","category":"page"},{"location":"man/formulations/#Specifying-inputs-in-NetworkStack","page":"Polar formulation","title":"Specifying inputs in NetworkStack","text":"","category":"section"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"Our three inputs are (v theta p_g) in mathbbR^2n_b + n_g (voltage magnitude, voltage angle, power generations). The basis psi is considered as an intermediate expression.","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"We store all inputs in a NetworkStack structure:","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"struct NetworkStack{VT} <: AbstractStack\n input::VT\n vmag::VT # voltage magnitudes (view)\n vang::VT # voltage angles (view)\n pgen::VT # active power generations (view)\n ψ::VT # nonlinear basis ψ(vmag, vang)\nend","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"All the inputs are specified in the vector input. The three vectors vmag, vang and pgen are views porting on input, and are defined mostly for convenience. By convention the vector input is ordered as [vmag; vang; pgen]:","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"# Define dimension of the problem\njulia> nbus, ngen, nlines = 3, 2, 4;\n\njulia> stack = ExaPF.NetworkStack(nbus, ngen, nlines, 0, Vector{Float64}, Vector{Float64})\n8-elements NetworkStack{Vector{Float64}}\n\njulia> stack.input\n8-element Vector{Float64}:\n 0.0\n 0.0\n 0.0\n 0.0\n 0.0\n 0.0\n 0.0\n 0.0\n\njulia> stack.vmag .= 1;\n\njulia> stack.vang .= 2;\n\njulia> stack.pgen .= 3;\n\njulia> stack.input\n8-element Vector{Float64}:\n 1.0\n 1.0\n 1.0\n 2.0\n 2.0\n 2.0\n 3.0\n 3.0","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"The basis vector ψ is an intermediate expression, whose values depend on the inputs.","category":"page"},{"location":"man/formulations/#Defining-a-state-and-a-control","page":"Polar formulation","title":"Defining a state and a control","text":"","category":"section"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"In the reduced space method, we have to split the variables in a state x and a control u. By default, we define","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":" x = (theta_pv theta_pq v_pq) quad\n x = (v_ref v_pv p_ggenpv) ","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"and the control, and was not flexible. In the new implementation, we define the state and the control as two mappings porting on the vector stack.input (which itself stores all the inputs in the problem):","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"julia> nbus, ngen, nlines = 4, 3, 4;\n\njulia> stack = ExaPF.NetworkStack(nbus, ngen, nlines, 0, Vector{Float64}, Vector{Float64});\n\njulia> stack.input .= 1:length(stack.input); # index array input\n\njulia> ref, pv, pq, genpv = [1], [2], [3, 4], [2, 3];\n\njulia> mapx = [nbus .+ pv; nbus .+ pq; pq];\n\njulia> mapu = [ref; pv; 2*nbus .+ genpv];\n\njulia> x = @view stack.input[mapx]\n5-element view(::Vector{Float64}, [6, 7, 8, 3, 4]) with eltype Float64:\n 6.0\n 7.0\n 8.0\n 3.0\n 4.0\n\njulia> u = @view stack.input[mapu]\n4-element view(::Vector{Float64}, [1, 2, 10, 11]) with eltype Float64:\n 1.0\n 2.0\n 10.0\n 11.0","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"By doing so, the values of the state and the control are directly stored inside the NetworkStack structure, avoiding to duplicate values in the memory.","category":"page"},{"location":"man/formulations/#How-to-manipulate-the-expressions?","page":"Polar formulation","title":"How to manipulate the expressions?","text":"","category":"section"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"ExaPF implements the different functions required to implement the optimal power flow problem with the polar formulation:","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"PowerFlowBalance: power flow balance equations\nPowerGenerationBounds: bounds on the power generation\nLineFlows: line flow constraints\nCostFunction: quadratic cost function","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"Each function follows the LTMR2020 model and depends on the basis function psi(v theta), here implemented in the PolarBasis function.","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"We demonstrate how to use the different functions on the case9 instance. The procedure remains the same for all power network.","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"julia> polar = ExaPF.load_polar(\"case9.m\");\n\njulia> stack = ExaPF.NetworkStack(polar);\n","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"note: Note\nAll the code presented below is agnostic with regards to the specific device (CPU, CUDABackend...) we are using. By default, ExaPF computes the expressions on the CPU. Deporting the computation on a CUDABackend simply translates to instantiate the PolarForm structure on the GPU: polar = PolarForm(\"case9.m\", CUDABackend()).","category":"page"},{"location":"man/formulations/#Interface","page":"Polar formulation","title":"Interface","text":"","category":"section"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"All functions are following AutoDiff.AbstractExpression's interface. The structure of the network is specified by the PolarForm we pass as an argument in the constructor. For instance, we build a new PolarBasis expression associated to case9 directly as","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"julia> basis = ExaPF.PolarBasis(polar)\nPolarBasis (AbstractExpression)\n","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"Each expression as a given dimension, given by","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"julia> length(basis)\n27\n","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"In ExaPF, the inputs and the parameters are stored inside a NetworkStack structure. Evaluating the basis psi naturally translates to","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"julia> basis(stack)\n27-element Vector{Float64}:\n 1.0\n 1.0\n 1.0\n 1.0\n 1.0\n 1.0\n 1.0\n 1.0\n 1.0\n 0.0\n ⋮\n 1.0\n 1.0\n 1.0\n 1.0\n 1.0\n 1.0\n 1.0\n 1.0\n 1.0","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"This function call allocates a vector psi with 27 elements and evaluates the basis associated to the LTMR2020 model. To avoid unnecessary allocation, one can preallocate the vector psi:","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"julia> psi = zeros(length(basis)) ;\n\njulia> basis(psi, stack);\n","category":"page"},{"location":"man/formulations/#Compose-expressions-together","page":"Polar formulation","title":"Compose expressions together","text":"","category":"section"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"In the LTMR2020 model, the polar basis psi(v theta) depends only on the voltage magnitudes and the voltage angles. However, this is not the case for the other functions (power flow balance, line flows, ...), which all depends on the basis psi(v theta).","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"In ExaPF, one has to build manually the vectorized expression tree associated to the power flow model. Luckily, evaluating the LTMR2020 simply amounts to compose functions together with the polar basis psi(v theta). ExaPF overloads the function ∘ to compose functions with a PolarBasis instance. The power flow balance can be evaluated as","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"julia> pflow = ExaPF.PowerFlowBalance(polar) ∘ basis;\n","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"which returns a ComposedExpressions structure.","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"The function pflow follows the same API, as any regular AutoDiff.AbstractExpression.","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"julia> n_balance = length(pflow)\n14\n\njulia> round.(pflow(stack); digits=6) # evaluate the power flow balance\n14-element Vector{Float64}:\n -1.63\n -0.85\n 0.0\n 0.9\n 0.0\n 1.0\n 0.0\n 1.25\n -0.167\n 0.042\n -0.2835\n 0.171\n -0.2275\n 0.259\n","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"When we evaluate a ComposedExpressions, ExaPF first computes the basis psi(v theta) inside stack.psi, and then ExaPF uses the values in stack.psi to evaluate the final result.","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"The procedure remains the same if one wants to evaluate the LineFlows or the PowerGenerationBounds. For instance, evaluating the line flows amounts to","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"julia> line_flows = ExaPF.LineFlows(polar) ∘ basis;\n\njulia> line_flows(stack)\n18-element Vector{Float64}:\n 0.0\n 0.006241000000000099\n 0.0320410000000001\n 0.0\n 0.010920249999999961\n 0.005550250000000068\n 0.0\n 0.02340899999999987\n 0.007743999999999858\n 0.0\n 0.006241000000000099\n 0.0320410000000001\n 0.0\n 0.010920249999999961\n 0.005550250000000068\n 0.0\n 0.02340899999999987\n 0.007743999999999858\n","category":"page"},{"location":"quickstart/#Quick-Start","page":"Quick start","title":"Quick Start","text":"","category":"section"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"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.","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"We start by importing CUDA and KernelAbstractions:","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"using CUDA\nusing KernelAbstractions","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"Then, we load ExaPF and its submodules with","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"using ExaPF\nimport ExaPF: AutoDiff\nconst PS = ExaPF.PowerSystem\nconst LS = ExaPF.LinearSolvers","category":"page"},{"location":"quickstart/#Short-version","page":"Quick start","title":"Short version","text":"","category":"section"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"ExaPF loads instances from the pglib-opf benchmark. ExaPF contains an artifact defined in Artifacts.toml that is built from the ExaData repository containing Exascale Computing Project relevant test cases. You may set a data file using","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"datafile = joinpath(artifact\"ExaData\", \"ExaData\", \"case1354.m\")","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"using LazyArtifacts\nusing ExaPF\nusing CUDA\nusing KernelAbstractions\nusing ExaPF\nimport ExaPF: AutoDiff\nconst PS = ExaPF.PowerSystem\nconst LS = ExaPF.LinearSolvers\nartifact_toml = joinpath(@__DIR__, \"..\", \"..\", \"Artifacts.toml\")\nexadata_hash = artifact_hash(\"ExaData\", artifact_toml)\ndatafile = joinpath(artifact_path(exadata_hash), \"ExaData\", \"case1354.m\")","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"The powerflow equations can be solved in three lines of code, as","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"polar = ExaPF.PolarForm(datafile, CPU()) # Load data\nstack = ExaPF.NetworkStack(polar) # Load variables\nconvergence = run_pf(polar, stack; verbose=1)\n","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"Implicitly, ExaPF has just proceed to the following operations:","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"instantiate automatically a starting point x_0 from the variables stored in stack\ninstantiate the Jacobian of the powerflow equations using AutoDiff.\nsolve the powerflow equations iteratively, using a Newton-Raphson algorithm.","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"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.","category":"page"},{"location":"quickstart/#Detailed-version","page":"Quick start","title":"Detailed version","text":"","category":"section"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"In what follows, we detail step by step the detailed procedure to solve the powerflow equations.","category":"page"},{"location":"quickstart/#How-to-load-a-MATPOWER-instance-as-a-PowerNetwork-object?","page":"Quick start","title":"How to load a MATPOWER instance as a PowerNetwork object?","text":"","category":"section"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"We start by importing a MATPOWER instance to a ExaPF.PowerSystem.PowerNetwork object:","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"pf = PS.PowerNetwork(datafile)","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"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","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"nbus = PS.get(pf, PS.NumberOfBuses())\npv_indexes = pf.pv;","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"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.","category":"page"},{"location":"quickstart/#How-to-solve-the-powerflow-equations?","page":"Quick start","title":"How to solve the powerflow equations?","text":"","category":"section"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"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 bmv = ve^j theta, each bus i=1 cdots N_B must satisfy the power balance equations:","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"beginaligned\n p_i = v_i sum_j^n v_j (g_ijcos(theta_i - theta_j) + b_ijsin(theta_i - theta_j)) \n q_i = v_i sum_j^n v_j (g_ijsin(theta_i - theta_j) - b_ijcos(theta_i - theta_j)) \nendaligned","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"The powerflow equations rewrite in the abstract mathematical formalism:","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"g(x u) = 0","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"For a given control u, solving the powerflow equations resumes to find a state x(u) such that g(x(u) u) = 0.","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"To this goal, ExaPF.jl implements a Newton-Raphson algorithm that allows to solve the powerflow equations in a few lines of code. We first instantiate a PolarForm object to adopt a polar formulation as a model:","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"polar = ExaPF.PolarForm(pf, CPU())","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"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().","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"The Newton-Raphson solves the equation g(x u) = 0 in an iterative fashion. The algorithm solves at each step the linear equation:","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":" x_k+1 = - (nabla_x g_k)^-1 g(x_k u)","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"Hence, the algorithm requires the following elements:","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"an initial variable x_0\na function to solve efficiently the linear system (nabla_x g_k) x_k+1 = g(x_k u)\na function to evaluate the Jacobian nabla_x g_k","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"The variable x is instantiated as:","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"stack = ExaPF.NetworkStack(polar)","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"The function g is implemented using ExaPF's custom modeler:","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"basis = ExaPF.PolarBasis(polar)\npowerflow = ExaPF.PowerFlowBalance(polar) ∘ basis","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"The Jacobian nabla_x g is evaluated automatically using forward-mode AutoDiff:","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"mapx = ExaPF.mapping(polar, State());\njx = ExaPF.Jacobian(polar, powerflow, mapx)","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"The (direct) linear solver can be instantiated directly as","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"linear_solver = LS.DirectSolver(jx.J);\n","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"Let's explain further these three objects.","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"stack is a AbstractStack storing all the variables attached to the formulation polar::PolarForm.\njx is a Jacobian structure which allows the solver to compute efficiently the Jacobian of the powerflow equations nabla_x g using AutoDiff.\nlinear_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.","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"In the AutoDiff Jacobian jx, the evaluation of the Jacobian J is stored in jx.J:","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"jac = jx.J","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"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.","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"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:","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"pf_algo = NewtonRaphson(; verbose=1, tol=1e-10)","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"Then, we can solve the powerflow equations simply with","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"convergence = ExaPF.nlsolve!(pf_algo, jx, stack; linear_solver=linear_solver)","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"Here, the algorithm solves the powerflow equations in 5 iterations. The algorithm modifies the values of stack inplace, to avoid any unnecessary memory allocations.","category":"page"},{"location":"quickstart/#How-to-deport-the-computation-on-the-GPU?","page":"Quick start","title":"How to deport the computation on the GPU?","text":"","category":"section"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"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:","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"polar_gpu = ExaPF.PolarForm(pf, CUDABackend())\n","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"polar_gpu will load all the structures it needs on the GPU, to avoid unnecessary movements between the host and the device. We can load the other structures directly on the GPU with:","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"stack_gpu = ExaPF.NetworkStack(polar_gpu)\n\nbasis_gpu = ExaPF.PolarBasis(polar_gpu)\npflow_gpu = ExaPF.PowerFlowBalance(polar_gpu) ∘ basis_gpu\njx_gpu = ExaPF.Jacobian(polar_gpu, pflow_gpu, mapx)\n\nlinear_solver = LS.DirectSolver(jx_gpu.J)","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"Then, solving the powerflow equations on the GPU directly translates as","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"convergence = ExaPF.nlsolve!(pf_algo, jx_gpu, stack_gpu; linear_solver=linear_solver)","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"Note that we get exactly the same iterations as when we solve the power flow equations on the CPU.","category":"page"},{"location":"quickstart/#How-to-solve-the-linear-system-with-BICGSTAB?","page":"Quick start","title":"How to solve the linear system with BICGSTAB?","text":"","category":"section"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"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).","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"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","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"npartitions = 8;\njac_gpu = jx_gpu.J;\nprecond = BlockJacobiPreconditioner(jac_gpu, npartitions, CUDABackend());","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"You can attach the preconditioner to an BICGSTAB algorithm simply as","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"linear_solver = ExaPF.KrylovBICGSTAB(jac_gpu; P=precond);\n","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"(this will use the BICGSTAB algorithm implemented in Krylov.jl).","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"We need to update accordingly the tolerance of the Newton-Raphson algorithm (the iterative solver is less accurate than the direct solver):","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"pf_algo = NewtonRaphson(; verbose=1, tol=1e-7)","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"We reset the variables to their initial values:","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"ExaPF.init!(polar_gpu, stack_gpu)","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"Then, solving the power flow with the iterative solvers directly translates to one call to nlsolve!:","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"convergence = ExaPF.nlsolve!(pf_algo, jx_gpu, stack_gpu; linear_solver=linear_solver)","category":"page"},{"location":"tutorials/batch_evaluation/","page":"Power flow: batch evaluation","title":"Power flow: batch evaluation","text":"CurrentModule = ExaPF\nDocTestSetup = quote\n using ExaPF\n const LS = ExaPF.LinearSolvers\n const AD = ExaPF.AD\nend\nDocTestFilters = [r\"ExaPF\"]","category":"page"},{"location":"tutorials/batch_evaluation/","page":"Power flow: batch evaluation","title":"Power flow: batch evaluation","text":"using ExaPF\nusing KLU\nusing LinearAlgebra\nconst LS = ExaPF.LinearSolvers\nconst PS = ExaPF.PowerSystem\n\npolar = ExaPF.load_polar(\"case9.m\")","category":"page"},{"location":"tutorials/batch_evaluation/#Batch-power-flow","page":"Power flow: batch evaluation","title":"Batch power flow","text":"","category":"section"},{"location":"tutorials/batch_evaluation/","page":"Power flow: batch evaluation","title":"Power flow: batch evaluation","text":"ExaPF provides a way to evaluate the expressions by blocks, opening the way to introduce more parallelism in the code.","category":"page"},{"location":"tutorials/batch_evaluation/#BlockPolarForm","page":"Power flow: batch evaluation","title":"BlockPolarForm","text":"","category":"section"},{"location":"tutorials/batch_evaluation/","page":"Power flow: batch evaluation","title":"Power flow: batch evaluation","text":"We recall that a given NetworkStack stack stores the different variables and parameters (power generations, voltages, loads) required to evaluate the power flow model.","category":"page"},{"location":"tutorials/batch_evaluation/","page":"Power flow: batch evaluation","title":"Power flow: batch evaluation","text":"stack = ExaPF.NetworkStack(polar);\n","category":"page"},{"location":"tutorials/batch_evaluation/","page":"Power flow: batch evaluation","title":"Power flow: batch evaluation","text":"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","category":"page"},{"location":"tutorials/batch_evaluation/","page":"Power flow: batch evaluation","title":"Power flow: batch evaluation","text":"nbus = ExaPF.get(polar, PS.NumberOfBuses());\npd = stack.params[1:nbus]\nqd = stack.params[nbus+1:2*nbus]\n","category":"page"},{"location":"tutorials/batch_evaluation/","page":"Power flow: batch evaluation","title":"Power flow: batch evaluation","text":"By default, a NetworkStack stores one set of loads p_0.","category":"page"},{"location":"tutorials/batch_evaluation/","page":"Power flow: batch evaluation","title":"Power flow: batch evaluation","text":"Suppose now we want to evaluate the model associated with the polar formulation for N different set of parameters (=scenarios) p_1 cdots p_N. ExaPF allows to streamline the polar formulation with a BlockPolarForm structure:","category":"page"},{"location":"tutorials/batch_evaluation/","page":"Power flow: batch evaluation","title":"Power flow: batch evaluation","text":"nscen = 10;\nblk_polar = ExaPF.BlockPolarForm(polar, nscen)\n","category":"page"},{"location":"tutorials/batch_evaluation/","page":"Power flow: batch evaluation","title":"Power flow: batch evaluation","text":"Then, ExaPF can also instantiate a NetworkStack object, with the memory required to store the variables of the different scenarios:","category":"page"},{"location":"tutorials/batch_evaluation/","page":"Power flow: batch evaluation","title":"Power flow: batch evaluation","text":"blk_stack = ExaPF.NetworkStack(blk_polar)\n","category":"page"},{"location":"tutorials/batch_evaluation/","page":"Power flow: batch evaluation","title":"Power flow: batch evaluation","text":"We can pass the scenarios manually using the function set_params!:","category":"page"},{"location":"tutorials/batch_evaluation/","page":"Power flow: batch evaluation","title":"Power flow: batch evaluation","text":"\nploads = rand(nbus, nscen);\nqloads = rand(nbus, nscen);\nExaPF.set_params!(blk_stack, ploads, qloads)\n","category":"page"},{"location":"tutorials/batch_evaluation/","page":"Power flow: batch evaluation","title":"Power flow: batch evaluation","text":"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):","category":"page"},{"location":"tutorials/batch_evaluation/","page":"Power flow: batch evaluation","title":"Power flow: batch evaluation","text":"reshape(blk_stack.vmag, nbus, nscen)","category":"page"},{"location":"tutorials/batch_evaluation/","page":"Power flow: batch evaluation","title":"Power flow: batch evaluation","text":"Only the parameters are varying according to the scenarios we passed as input in the constructor:","category":"page"},{"location":"tutorials/batch_evaluation/","page":"Power flow: batch evaluation","title":"Power flow: batch evaluation","text":"reshape(blk_stack.pload, nbus, nscen)","category":"page"},{"location":"tutorials/batch_evaluation/#Evaluate-expressions-in-block","page":"Power flow: batch evaluation","title":"Evaluate expressions in block","text":"","category":"section"},{"location":"tutorials/batch_evaluation/","page":"Power flow: batch evaluation","title":"Power flow: batch evaluation","text":"ExaPF takes advantage of the block structure when using a BlockPolarForm.","category":"page"},{"location":"tutorials/batch_evaluation/","page":"Power flow: batch evaluation","title":"Power flow: batch evaluation","text":"As an example, suppose we want to evaluate the power flow balances in block form with a PowerFlowBalance expression:","category":"page"},{"location":"tutorials/batch_evaluation/","page":"Power flow: batch evaluation","title":"Power flow: batch evaluation","text":"powerflow = ExaPF.PowerFlowBalance(blk_polar) ∘ ExaPF.PolarBasis(blk_polar);\n","category":"page"},{"location":"tutorials/batch_evaluation/","page":"Power flow: batch evaluation","title":"Power flow: batch evaluation","text":"A block evaluation takes as input the NetworkStack blk_stack structure:","category":"page"},{"location":"tutorials/batch_evaluation/","page":"Power flow: batch evaluation","title":"Power flow: batch evaluation","text":"m = div(length(powerflow), nscen);\nblk_output = zeros(m * nscen);\npowerflow(blk_output, blk_stack); # inplace evaluation\nreshape(blk_output, m, nscen)\n","category":"page"},{"location":"tutorials/batch_evaluation/","page":"Power flow: batch evaluation","title":"Power flow: batch evaluation","text":"We get N different results for the power flow balance equations, depending on which scenario we are on.","category":"page"},{"location":"tutorials/batch_evaluation/#Solve-power-flow-in-block-on-the-CPU","page":"Power flow: batch evaluation","title":"Solve power flow in block on the CPU","text":"","category":"section"},{"location":"tutorials/batch_evaluation/","page":"Power flow: batch evaluation","title":"Power flow: batch evaluation","text":"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:","category":"page"},{"location":"tutorials/batch_evaluation/","page":"Power flow: batch evaluation","title":"Power flow: batch evaluation","text":"blk_jx = ExaPF.ArrowheadJacobian(blk_polar, powerflow, State());\nblk_jx.J","category":"page"},{"location":"tutorials/batch_evaluation/","page":"Power flow: batch evaluation","title":"Power flow: batch evaluation","text":"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:","category":"page"},{"location":"tutorials/batch_evaluation/","page":"Power flow: batch evaluation","title":"Power flow: batch evaluation","text":"ExaPF.set_params!(blk_jx, blk_stack);\nExaPF.jacobian!(blk_jx, blk_stack);\n","category":"page"},{"location":"tutorials/batch_evaluation/","page":"Power flow: batch evaluation","title":"Power flow: batch evaluation","text":"As soon as the blk_jx initialized, we can solve the power flow equations in block as","category":"page"},{"location":"tutorials/batch_evaluation/","page":"Power flow: batch evaluation","title":"Power flow: batch evaluation","text":"conv = ExaPF.nlsolve!(\n NewtonRaphson(verbose=2),\n blk_jx,\n blk_stack;\n)","category":"page"},{"location":"tutorials/batch_evaluation/","page":"Power flow: batch evaluation","title":"Power flow: batch evaluation","text":"At the solution, we get different values for the voltage magnitudes at the PQ nodes:","category":"page"},{"location":"tutorials/batch_evaluation/","page":"Power flow: batch evaluation","title":"Power flow: batch evaluation","text":"reshape(blk_stack.vmag, nbus, nscen)","category":"page"},{"location":"tutorials/batch_evaluation/#Solve-power-flow-in-batch-on-the-GPU","page":"Power flow: batch evaluation","title":"Solve power flow in batch on the GPU","text":"","category":"section"},{"location":"tutorials/batch_evaluation/","page":"Power flow: batch evaluation","title":"Power flow: batch evaluation","text":"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:","category":"page"},{"location":"tutorials/batch_evaluation/","page":"Power flow: batch evaluation","title":"Power flow: batch evaluation","text":"using CUDA, CUSOLVERRF\npolar_gpu = ExaPF.load_polar(\"case9.m\", CUDABackend());\nblk_polar_gpu = ExaPF.BlockPolarForm(polar_gpu, nscen); # load model on GPU\nblk_stack_gpu = ExaPF.NetworkStack(blk_polar_gpu);\nExaPF.set_params!(blk_stack_gpu, ploads, qloads);\npowerflow_gpu = ExaPF.PowerFlowBalance(blk_polar_gpu) ∘ ExaPF.PolarBasis(blk_polar_gpu);\nblk_jx_gpu = ExaPF.ArrowheadJacobian(blk_polar_gpu, powerflow_gpu, State());\nExaPF.set_params!(blk_jx_gpu, blk_stack_gpu);\nExaPF.jacobian!(blk_jx_gpu, blk_stack_gpu);\nrf_fac = CUSOLVERRF.RFLU(blk_jx_gpu.J)\nrf_solver = LS.DirectSolver(rf_fac)\nconv = ExaPF.nlsolve!(\n NewtonRaphson(verbose=2),\n blk_jx_gpu,\n blk_stack_gpu;\n linear_solver=rf_solver,\n)\n","category":"page"},{"location":"lib/powersystem/","page":"PowerSystem","title":"PowerSystem","text":"CurrentModule = ExaPF.PowerSystem","category":"page"},{"location":"lib/powersystem/#PowerSystem","page":"PowerSystem","title":"PowerSystem","text":"","category":"section"},{"location":"lib/powersystem/#Description","page":"PowerSystem","title":"Description","text":"","category":"section"},{"location":"lib/powersystem/","page":"PowerSystem","title":"PowerSystem","text":"AbstractPowerSystem\nPowerNetwork\nload_case","category":"page"},{"location":"lib/powersystem/#ExaPF.PowerSystem.AbstractPowerSystem","page":"PowerSystem","title":"ExaPF.PowerSystem.AbstractPowerSystem","text":"AbstractPowerSystem\n\nFirst layer of the package. Store the topology of a given transmission network, including:\n\nthe power injection at each bus ;\nthe admittance matrix ;\nthe default voltage at each bus.\n\nData are imported either from a matpower file, or a PSSE file.\n\n\n\n\n\n","category":"type"},{"location":"lib/powersystem/#ExaPF.PowerSystem.PowerNetwork","page":"PowerSystem","title":"ExaPF.PowerSystem.PowerNetwork","text":"PowerNetwork <: AbstractPowerSystem\n\nThis structure contains constant parameters that define the topology and physics of the power network.\n\nThe 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.\n\nNote\n\nThe object PowerNetwork is created in the host memory. Use a AbstractFormulation object to move data to the target device.\n\n\n\n\n\n","category":"type"},{"location":"lib/powersystem/#ExaPF.PowerSystem.load_case","page":"PowerSystem","title":"ExaPF.PowerSystem.load_case","text":"load_case(case_name, lib::PowerNetworkLibrary=EXADATA)\n\nConvenient 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.\n\nExamples\n\njulia> net = PS.load_case(\"case118\") # default is MATPOWER\nPowerNetwork object with:\n Buses: 118 (Slack: 1. PV: 53. PQ: 64)\n Generators: 54.\n\njulia> net = PS.load_case(\"case1354_pegase\", PS.PGLIB)\nPowerNetwork object with:\n Buses: 1354 (Slack: 1. PV: 259. PQ: 1094)\n Generators: 260.\n\n\n\n\n\n","category":"function"},{"location":"lib/powersystem/#API-Reference","page":"PowerSystem","title":"API Reference","text":"","category":"section"},{"location":"lib/powersystem/#Network-elements","page":"PowerSystem","title":"Network elements","text":"","category":"section"},{"location":"lib/powersystem/","page":"PowerSystem","title":"PowerSystem","text":"AbstractNetworkElement","category":"page"},{"location":"lib/powersystem/#ExaPF.PowerSystem.AbstractNetworkElement","page":"PowerSystem","title":"ExaPF.PowerSystem.AbstractNetworkElement","text":"AbstractNetworkElement\n\nAbstraction for all physical elements being parts of a AbstractPowerSystem. Elements are divided in\n\ntransmission lines (Lines)\nbuses (Buses)\ngenerators (Generators)\n\n\n\n\n\n","category":"type"},{"location":"lib/powersystem/","page":"PowerSystem","title":"PowerSystem","text":"List of elements:","category":"page"},{"location":"lib/powersystem/","page":"PowerSystem","title":"PowerSystem","text":"Buses\nLines\nGenerators","category":"page"},{"location":"lib/powersystem/#ExaPF.PowerSystem.Buses","page":"PowerSystem","title":"ExaPF.PowerSystem.Buses","text":"Buses <: AbstractNetworkElement\n\nBuses of a transmission network.\n\n\n\n\n\n","category":"type"},{"location":"lib/powersystem/#ExaPF.PowerSystem.Lines","page":"PowerSystem","title":"ExaPF.PowerSystem.Lines","text":"Lines <: AbstractNetworkElement\n\nLines of a transmission network.\n\n\n\n\n\n","category":"type"},{"location":"lib/powersystem/#ExaPF.PowerSystem.Generators","page":"PowerSystem","title":"ExaPF.PowerSystem.Generators","text":"Generators <: AbstractElement\n\nGenerators in a transmission network\n\n\n\n\n\n","category":"type"},{"location":"lib/powersystem/#Network-attributes","page":"PowerSystem","title":"Network attributes","text":"","category":"section"},{"location":"lib/powersystem/","page":"PowerSystem","title":"PowerSystem","text":"AbstractNetworkAttribute","category":"page"},{"location":"lib/powersystem/#ExaPF.PowerSystem.AbstractNetworkAttribute","page":"PowerSystem","title":"ExaPF.PowerSystem.AbstractNetworkAttribute","text":"AbstractNetworkAttribute\n\nAttribute of a AbstractPowerSystem.\n\n\n\n\n\n","category":"type"},{"location":"lib/powersystem/","page":"PowerSystem","title":"PowerSystem","text":"List of attributes:","category":"page"},{"location":"lib/powersystem/","page":"PowerSystem","title":"PowerSystem","text":"NumberOfBuses\nNumberOfLines\nNumberOfGenerators\nNumberOfPVBuses\nNumberOfPQBuses\nNumberOfSlackBuses\nBaseMVA\nBusAdmittanceMatrix","category":"page"},{"location":"lib/powersystem/#ExaPF.PowerSystem.NumberOfBuses","page":"PowerSystem","title":"ExaPF.PowerSystem.NumberOfBuses","text":"NumberOfBuses <: AbstractNetworkAttribute\n\nNumber of buses in a AbstractPowerSystem.\n\n\n\n\n\n","category":"type"},{"location":"lib/powersystem/#ExaPF.PowerSystem.NumberOfLines","page":"PowerSystem","title":"ExaPF.PowerSystem.NumberOfLines","text":"NumberOfLines <: AbstractNetworkAttribute\n\nNumber of lines in a AbstractPowerSystem.\n\n\n\n\n\n","category":"type"},{"location":"lib/powersystem/#ExaPF.PowerSystem.NumberOfGenerators","page":"PowerSystem","title":"ExaPF.PowerSystem.NumberOfGenerators","text":"NumberOfGenerators <: AbstractNetworkAttribute\n\nNumber of generators in a AbstractPowerSystem.\n\n\n\n\n\n","category":"type"},{"location":"lib/powersystem/#ExaPF.PowerSystem.NumberOfPVBuses","page":"PowerSystem","title":"ExaPF.PowerSystem.NumberOfPVBuses","text":"NumberOfPVBuses <: AbstractNetworkAttribute\n\nNumber of PV buses in a AbstractPowerSystem.\n\n\n\n\n\n","category":"type"},{"location":"lib/powersystem/#ExaPF.PowerSystem.NumberOfPQBuses","page":"PowerSystem","title":"ExaPF.PowerSystem.NumberOfPQBuses","text":"NumberOfPQBuses <: AbstractNetworkAttribute\n\nNumber of PQ buses in a AbstractPowerSystem.\n\n\n\n\n\n","category":"type"},{"location":"lib/powersystem/#ExaPF.PowerSystem.NumberOfSlackBuses","page":"PowerSystem","title":"ExaPF.PowerSystem.NumberOfSlackBuses","text":"NumberOfSlackBuses <: AbstractNetworkAttribute\n\nNumber of slack buses in a AbstractPowerSystem.\n\n\n\n\n\n","category":"type"},{"location":"lib/powersystem/#ExaPF.PowerSystem.BaseMVA","page":"PowerSystem","title":"ExaPF.PowerSystem.BaseMVA","text":"BaseMVA <: AbstractNetworkAttribute\n\nBase MVA of the network.\n\n\n\n\n\n","category":"type"},{"location":"lib/powersystem/#ExaPF.PowerSystem.BusAdmittanceMatrix","page":"PowerSystem","title":"ExaPF.PowerSystem.BusAdmittanceMatrix","text":"BusAdmittanceMatrix <: AbstractNetworkAttribute\n\nBus admittance matrix associated with the topology of the network.\n\n\n\n\n\n","category":"type"},{"location":"lib/powersystem/","page":"PowerSystem","title":"PowerSystem","text":"Query the indexing of the different elements in a given network:","category":"page"},{"location":"lib/powersystem/","page":"PowerSystem","title":"PowerSystem","text":"PVIndexes\nPQIndexes\nSlackIndexes\nGeneratorIndexes\n","category":"page"},{"location":"lib/powersystem/#ExaPF.PowerSystem.PVIndexes","page":"PowerSystem","title":"ExaPF.PowerSystem.PVIndexes","text":"PVIndexes <: AbstractIndexing\n\nIndexes of the PV buses in a AbstractPowerSystem.\n\n\n\n\n\n","category":"type"},{"location":"lib/powersystem/#ExaPF.PowerSystem.PQIndexes","page":"PowerSystem","title":"ExaPF.PowerSystem.PQIndexes","text":"PQIndexes <: AbstractIndexing\n\nIndexes of the PQ buses in a AbstractPowerSystem.\n\n\n\n\n\n","category":"type"},{"location":"lib/powersystem/#ExaPF.PowerSystem.SlackIndexes","page":"PowerSystem","title":"ExaPF.PowerSystem.SlackIndexes","text":"SlackIndexes <: AbstractIndexing\n\nIndexes of the slack buses in a AbstractPowerSystem.\n\n\n\n\n\n","category":"type"},{"location":"lib/powersystem/#ExaPF.PowerSystem.GeneratorIndexes","page":"PowerSystem","title":"ExaPF.PowerSystem.GeneratorIndexes","text":"GeneratorIndexes <: AbstractIndexing\n\nIndexes of the generators in a AbstractPowerSystem.\n\n\n\n\n\n","category":"type"},{"location":"lib/powersystem/#Network-values","page":"PowerSystem","title":"Network values","text":"","category":"section"},{"location":"lib/powersystem/","page":"PowerSystem","title":"PowerSystem","text":"AbstractNetworkValues","category":"page"},{"location":"lib/powersystem/#ExaPF.PowerSystem.AbstractNetworkValues","page":"PowerSystem","title":"ExaPF.PowerSystem.AbstractNetworkValues","text":"AbstractNetworkValues\n\nNumerical values attached to the different attributes of the network.\n\n\n\n\n\n","category":"type"},{"location":"lib/powersystem/","page":"PowerSystem","title":"PowerSystem","text":"List of values:","category":"page"},{"location":"lib/powersystem/","page":"PowerSystem","title":"PowerSystem","text":"VoltageMagnitude\nVoltageAngle\nActivePower\nReactivePower\nActiveLoad\nReactiveLoad\n","category":"page"},{"location":"lib/powersystem/#ExaPF.PowerSystem.VoltageMagnitude","page":"PowerSystem","title":"ExaPF.PowerSystem.VoltageMagnitude","text":"VoltageMagnitude <: AbstractNetworkValues\n\nMagnitude |v| of the voltage v = |v| exp(i θ).\n\n\n\n\n\n","category":"type"},{"location":"lib/powersystem/#ExaPF.PowerSystem.VoltageAngle","page":"PowerSystem","title":"ExaPF.PowerSystem.VoltageAngle","text":"VoltageAngle <: AbstractNetworkValues\n\nAngle θ of the voltage v = |v| exp(i θ).\n\n\n\n\n\n","category":"type"},{"location":"lib/powersystem/#ExaPF.PowerSystem.ActivePower","page":"PowerSystem","title":"ExaPF.PowerSystem.ActivePower","text":"ActivePower <: AbstractNetworkValues\n\nActive power P of the complex power S = P + iQ.\n\n\n\n\n\n","category":"type"},{"location":"lib/powersystem/#ExaPF.PowerSystem.ReactivePower","page":"PowerSystem","title":"ExaPF.PowerSystem.ReactivePower","text":"ReactivePower <: AbstractNetworkValues\n\nReactive power Q of the complex power S = P + iQ.\n\n\n\n\n\n","category":"type"},{"location":"lib/powersystem/#ExaPF.PowerSystem.ActiveLoad","page":"PowerSystem","title":"ExaPF.PowerSystem.ActiveLoad","text":"ActiveLoad <: AbstractNetworkValues\n\nActive load Pd at buses.\n\n\n\n\n\n","category":"type"},{"location":"lib/powersystem/#ExaPF.PowerSystem.ReactiveLoad","page":"PowerSystem","title":"ExaPF.PowerSystem.ReactiveLoad","text":"ReactiveLoad <: AbstractNetworkValues\n\nReactive load Qd at buses.\n\n\n\n\n\n","category":"type"},{"location":"lib/powersystem/","page":"PowerSystem","title":"PowerSystem","text":"Function to get the range of a given value:","category":"page"},{"location":"lib/powersystem/","page":"PowerSystem","title":"PowerSystem","text":"bounds","category":"page"},{"location":"lib/powersystem/#ExaPF.PowerSystem.bounds","page":"PowerSystem","title":"ExaPF.PowerSystem.bounds","text":"bounds(pf::AbstractPowerSystem, attr::AbstractNetworkAttribute, val::AbstractNetworkValues)\n\nReturn lower and upper bounds corresponding to the admissible values of the AbstractNetworkAttribute attr.\n\nExamples\n\np_min, p_max = bounds(pf, Generator(), ActivePower())\nv_min, v_max = bounds(pf, Buses(), VoltageMagnitude())\n\n\n\n\n\n\n","category":"function"},{"location":"#ExaPF","page":"Home","title":"ExaPF","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"ExaPF.jl is a package to solve the power flow problem on upcoming exascale architectures. The code has been designed to be:","category":"page"},{"location":"","page":"Home","title":"Home","text":"Portable: Targeting exascale architectures implies a focus on graphics processing units (GPUs) as these systems lack substantial computational performance through classical CPUs.\nDifferentiable: All the expressions implemented in ExaPF are fully compatible with ForwardDiff.jl, and routines are provided to extract first- and second-order derivatives to solve efficiently power flow and optimal power flow problems.","category":"page"},{"location":"","page":"Home","title":"Home","text":"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 provide extensions that leverage the packages CUDA.jl, [AMDGPU.jl]((https://github.com/JuliaGPU/AMDGPU.jl), and KernelAbstractions.jl to make ExaPF portable across GPU architectures.","category":"page"},{"location":"#Table-of-contents","page":"Home","title":"Table of contents","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"Pages = [\n \"quickstart.md\",\n]\nDepth=1","category":"page"},{"location":"#Manual","page":"Home","title":"Manual","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"Pages = [\n \"man/formulations.md\",\n \"man/powersystem.md\",\n \"man/autodiff.md\",\n \"man/linearsolver.md\",\n \"man/benchmark.md\",\n]\nDepth = 1","category":"page"},{"location":"#Library","page":"Home","title":"Library","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"Pages = [\n \"lib/formulations.md\",\n \"lib/powersystem.md\",\n \"lib/autodiff.md\",\n \"lib/linearsolver.md\",\n]\nDepth = 1","category":"page"},{"location":"#Artifact","page":"Home","title":"Artifact","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"Pages = [\n \"artifact.md\",\n]\nDepth = 1","category":"page"},{"location":"#Funding","page":"Home","title":"Funding","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"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.","category":"page"},{"location":"lib/linearsolver/","page":"Linear Solvers","title":"Linear Solvers","text":"CurrentModule = ExaPF.LinearSolvers","category":"page"},{"location":"lib/linearsolver/#Linear-solvers","page":"Linear Solvers","title":"Linear solvers","text":"","category":"section"},{"location":"lib/linearsolver/#Description","page":"Linear Solvers","title":"Description","text":"","category":"section"},{"location":"lib/linearsolver/","page":"Linear Solvers","title":"Linear Solvers","text":"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!.","category":"page"},{"location":"lib/linearsolver/","page":"Linear Solvers","title":"Linear Solvers","text":"ldiv!","category":"page"},{"location":"lib/linearsolver/#ExaPF.LinearSolvers.ldiv!","page":"Linear Solvers","title":"ExaPF.LinearSolvers.ldiv!","text":"ldiv!(solver, x, A, y)\nldiv!(solver, x, y)\n\nsolver::AbstractLinearSolver: linear solver to solve the system\nx::AbstractVector: Solution\nA::AbstractMatrix: Input matrix\ny::AbstractVector: RHS\n\nSolve 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.\n\n\n\n\n\n","category":"function"},{"location":"lib/linearsolver/#Direct-solvers","page":"Linear Solvers","title":"Direct solvers","text":"","category":"section"},{"location":"lib/linearsolver/","page":"Linear Solvers","title":"Linear Solvers","text":"ExaPF wraps UMFPACK (shipped with SuiteSparse.jl) on the CPU, and CUSPARSE on CUDA device.","category":"page"},{"location":"lib/linearsolver/","page":"Linear Solvers","title":"Linear Solvers","text":"DirectSolver","category":"page"},{"location":"lib/linearsolver/#ExaPF.LinearSolvers.DirectSolver","page":"Linear Solvers","title":"ExaPF.LinearSolvers.DirectSolver","text":"DirectSolver <: AbstractLinearSolver\n\nSolve linear system A x = y with direct linear algebra.\n\nOn the CPU, DirectSolver uses UMFPACK to solve the linear system\nOn CUDA GPU, DirectSolver redirects the resolution to the method CUSOLVER.csrlsvqr\n\n\n\n\n\n","category":"type"},{"location":"lib/linearsolver/#Iterative-solvers","page":"Linear Solvers","title":"Iterative solvers","text":"","category":"section"},{"location":"lib/linearsolver/","page":"Linear Solvers","title":"Linear Solvers","text":"KrylovBICGSTAB\nDQGMRES\nBICGSTAB\nEigenBICGSTAB","category":"page"},{"location":"lib/linearsolver/#ExaPF.LinearSolvers.KrylovBICGSTAB","page":"Linear Solvers","title":"ExaPF.LinearSolvers.KrylovBICGSTAB","text":"KrylovBICGSTAB <: AbstractIterativeLinearSolver\nKrylovBICGSTAB(precond; verbose=0, rtol=1e-10, atol=1e-10)\n\nWrap Krylov.jl's BICGSTAB algorithm to solve iteratively the linear system A x = y.\n\n\n\n\n\n","category":"type"},{"location":"lib/linearsolver/#ExaPF.LinearSolvers.DQGMRES","page":"Linear Solvers","title":"ExaPF.LinearSolvers.DQGMRES","text":"DQGMRES <: AbstractIterativeLinearSolver\nDQGMRES(precond; verbose=false, memory=4)\n\nWrap Krylov.jl's DQGMRES algorithm to solve iteratively the linear system A x = y.\n\n\n\n\n\n","category":"type"},{"location":"lib/linearsolver/#ExaPF.LinearSolvers.BICGSTAB","page":"Linear Solvers","title":"ExaPF.LinearSolvers.BICGSTAB","text":"BICGSTAB <: AbstractIterativeLinearSolver\nBICGSTAB(precond; maxiter=2_000, tol=1e-8, verbose=false)\n\nCustom BICGSTAB implementation to solve iteratively the linear system A x = y.\n\n\n\n\n\n","category":"type"},{"location":"lib/linearsolver/#ExaPF.LinearSolvers.EigenBICGSTAB","page":"Linear Solvers","title":"ExaPF.LinearSolvers.EigenBICGSTAB","text":"EigenBICGSTAB <: AbstractIterativeLinearSolver\nEigenBICGSTAB(precond; maxiter=2_000, tol=1e-8, verbose=false)\n\nJulia's port of Eigen's BICGSTAB to solve iteratively the linear system A x = y.\n\n\n\n\n\n","category":"type"},{"location":"lib/linearsolver/","page":"Linear Solvers","title":"Linear Solvers","text":"ExaPF.jl is shipped with a custom BICGSTAB implementation. However, we highly recommend to use KrylovBICGSTAB instead, which has proved to be more robust.","category":"page"},{"location":"lib/linearsolver/","page":"Linear Solvers","title":"Linear Solvers","text":"bicgstab\n","category":"page"},{"location":"lib/linearsolver/#ExaPF.LinearSolvers.bicgstab","page":"Linear Solvers","title":"ExaPF.LinearSolvers.bicgstab","text":"bicgstab(A, b, P, xi;\n tol=1e-8,\n maxiter=size(A, 1),\n verbose=false,\n maxtol=1e20)\n\nBiCGSTAB implementation according to\n\nVan 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.\n\n\n\n\n\n","category":"function"},{"location":"lib/linearsolver/","page":"Linear Solvers","title":"Linear Solvers","text":"Available linear solvers can be queried with","category":"page"},{"location":"lib/linearsolver/","page":"Linear Solvers","title":"Linear Solvers","text":"list_solvers\n","category":"page"},{"location":"lib/linearsolver/#ExaPF.LinearSolvers.list_solvers","page":"Linear Solvers","title":"ExaPF.LinearSolvers.list_solvers","text":"list_solvers(::KernelAbstractions.Device)\n\nList linear solvers available on current device.\n\n\n\n\n\n","category":"function"},{"location":"lib/linearsolver/","page":"Linear Solvers","title":"Linear Solvers","text":"A default solver is provided for each vendor backend.","category":"page"},{"location":"lib/linearsolver/","page":"Linear Solvers","title":"Linear Solvers","text":"default_linear_solver\n","category":"page"},{"location":"lib/linearsolver/#ExaPF.LinearSolvers.default_linear_solver","page":"Linear Solvers","title":"ExaPF.LinearSolvers.default_linear_solver","text":"default_linear_solver(A, ::KA.CPU)\n\nDefault linear solver on the CPU.\n\n\n\n\n\n","category":"function"},{"location":"lib/linearsolver/#Block-Krylov-solvers","page":"Linear Solvers","title":"Block-Krylov solvers","text":"","category":"section"},{"location":"lib/linearsolver/","page":"Linear Solvers","title":"Linear Solvers","text":"ExaPF.jl provides an implementation of block-GMRES to solve linear systems with multiple righ-hand sides.","category":"page"},{"location":"lib/linearsolver/","page":"Linear Solvers","title":"Linear Solvers","text":"BlockKrylovSolver\nBlockGmresSolver\nblock_gmres\nblock_gmres!","category":"page"},{"location":"lib/linearsolver/#ExaPF.LinearSolvers.BlockKrylovSolver","page":"Linear Solvers","title":"ExaPF.LinearSolvers.BlockKrylovSolver","text":"Abstract type for using block Krylov solvers in-place\n\n\n\n\n\n","category":"type"},{"location":"lib/linearsolver/#ExaPF.LinearSolvers.BlockGmresSolver","page":"Linear Solvers","title":"ExaPF.LinearSolvers.BlockGmresSolver","text":"Type for storing the vectors required by the in-place version of BLOCK-GMRES.\n\nThe outer constructors\n\nsolver = BlockGmresSolver(m, n, p, memory, SV, SM)\nsolver = BlockGmresSolver(A, B; memory=5)\n\nmay be used in order to create these vectors. memory is set to div(n,p) if the value given is larger than div(n,p).\n\n\n\n\n\n","category":"type"},{"location":"lib/linearsolver/#ExaPF.LinearSolvers.block_gmres","page":"Linear Solvers","title":"ExaPF.LinearSolvers.block_gmres","text":"(X, stats) = block_gmres(A, B; X0::AbstractMatrix{FC}; memory::Int=20, M=I, N=I,\n ldiv::Bool=false, restart::Bool=false, reorthogonalization::Bool=false,\n atol::T = √eps(T), rtol::T=√eps(T), itmax::Int=0,\n timemax::Float64=Inf, verbose::Int=0, history::Bool=false)\n\nT is an AbstractFloat such as Float32, Float64 or BigFloat. FC is T or Complex{T}.\n\n(X, stats) = block_gmres(A, B, X0::AbstractVector; kwargs...)\n\nGMRES can be warm-started from an initial guess X0 where kwargs are the same keyword arguments as above.\n\nSolve the linear system AX = B of size n with p right-hand sides using block-GMRES.\n\nInput arguments\n\nA: a linear operator that models a matrix of dimension n;\nB: a matrix of size n × p.\n\nOptional argument\n\nX0: a matrix of size n × p that represents an initial guess of the solution X.\n\nKeyword arguments\n\nmemory: if restart = true, the restarted version block-GMRES(k) is used with k = memory. If restart = false, the parameter memory should be used as a hint of the number of iterations to limit dynamic memory allocations. Additional storage will be allocated if the number of iterations exceeds memory;\nM: linear operator that models a nonsingular matrix of size n used for left preconditioning;\nN: linear operator that models a nonsingular matrix of size n used for right preconditioning;\nldiv: define whether the preconditioners use ldiv! or mul!;\nrestart: restart the method after memory iterations;\nreorthogonalization: reorthogonalize the new matrices of the block-Krylov basis against all previous matrix;\natol: absolute stopping tolerance based on the residual norm;\nrtol: relative stopping tolerance based on the residual norm;\nitmax: the maximum number of iterations. If itmax=0, the default number of iterations is set to 2 * div(n,p);\ntimemax: the time limit in seconds;\nverbose: additional details can be displayed if verbose mode is enabled (verbose > 0). Information will be displayed every verbose iterations;\nhistory: collect additional statistics on the run such as residual norms.\n\nOutput arguments\n\nx: a dense matrix of size n × p;\nstats: statistics collected on the run in a BlockGmresStats.\n\n\n\n\n\n","category":"function"},{"location":"lib/linearsolver/#ExaPF.LinearSolvers.block_gmres!","page":"Linear Solvers","title":"ExaPF.LinearSolvers.block_gmres!","text":"solver = block_gmres!(solver::BlockGmresSolver, B; kwargs...)\nsolver = block_gmres!(solver::BlockGmresSolver, B, X0; kwargs...)\n\nwhere kwargs are keyword arguments of block_gmres.\n\nSee BlockGmresSolver for more details about the solver.\n\n\n\n\n\n","category":"function"},{"location":"man/autodiff/","page":"AutoDiff","title":"AutoDiff","text":"DocTestSetup = quote\n using ExaPF\n const AD = ExaPF.AD\nend\nDocTestFilters = [r\"ExaPF\"]","category":"page"},{"location":"man/autodiff/#AutoDiff","page":"AutoDiff","title":"AutoDiff","text":"","category":"section"},{"location":"man/autodiff/#Overview","page":"AutoDiff","title":"Overview","text":"","category":"section"},{"location":"man/autodiff/","page":"AutoDiff","title":"AutoDiff","text":"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.","category":"page"},{"location":"man/autodiff/","page":"AutoDiff","title":"AutoDiff","text":" go = true\n while(go)\n dx .= jacobian(x)\\f(x)\n x .= x .- dx\n go = norm(f(x)) < tol ? true : false\n end","category":"page"},{"location":"man/autodiff/","page":"AutoDiff","title":"AutoDiff","text":"There are two modes of differentiation called forward/tangent or reverse/adjoint. The latter is known in machine learning as backpropagation. The forward mode generates Jacobian-vector product code tgt(x,d) = J(x) * d, while the adjoint mode generates code for the transposed Jacobian-vector product adj(x,y) = (J(x)'*y). We recommend the book Evaluating derivatives: principles and techniques of algorithmic differentiation by Griewank and Walther[1] for a more in-depth introduction to automatic differentiation. The computational complexity of both models favors the adjoint mode if the number of outputs of F is much smaller than the number of inputs size(x) >> size(F), like for example the loss functions in machine learning. However, in our case F is a multivariate vector function from mathbbR^n to mathbbR^n, where n is the number of buses.","category":"page"},{"location":"man/autodiff/","page":"AutoDiff","title":"AutoDiff","text":"(Image: Jacobian coloring \\label{fig:coloring})","category":"page"},{"location":"man/autodiff/","page":"AutoDiff","title":"AutoDiff","text":"To avoid a complexity of mathcalO(n) cdot cost(F) by letting the tangent mode run over all Cartesian basis vectors of mathbbR^n, we apply the technique of Jacobian coloring to compress the sparse Jacobian J. Running the tangent mode, it allows to compute columns of the Jacobian concurrently, by combining independent columns in one Jacobian-vector evaluation (see in figure above). For sparsity detection we rely on the greedy algorithm implemented by SparseDiffTools.jl.","category":"page"},{"location":"man/autodiff/","page":"AutoDiff","title":"AutoDiff","text":"Given the sparsity pattern, the forward model is applied through the package ForwardDiff.jl. Given the number of Jacobian colors c we can build our dual type t1s with c directions:","category":"page"},{"location":"man/autodiff/","page":"AutoDiff","title":"AutoDiff","text":"t1s{N} = ForwardDiff.Dual{Nothing,Float64, N} where N}","category":"page"},{"location":"man/autodiff/","page":"AutoDiff","title":"AutoDiff","text":"Note that a second-order type t2s can be created naturally by applying the same logic to t1s:","category":"page"},{"location":"man/autodiff/","page":"AutoDiff","title":"AutoDiff","text":"t2s{M,N} = ForwardDiff.Dual{Nothing,t1s{N}, M} where M, N}","category":"page"},{"location":"man/autodiff/","page":"AutoDiff","title":"AutoDiff","text":"Finally, this dual type can be ported to both vector types Vector and CuVector:","category":"page"},{"location":"man/autodiff/","page":"AutoDiff","title":"AutoDiff","text":"VT = Vector{Float64}\nVT = Vector{t1s{N}}}\nVT = CuVector{t1s{N}}}","category":"page"},{"location":"man/autodiff/","page":"AutoDiff","title":"AutoDiff","text":"Setting VT to either of the three types allows us to instantiate code that has been written using the broadcast operator .","category":"page"},{"location":"man/autodiff/","page":"AutoDiff","title":"AutoDiff","text":"x .= a .* b","category":"page"},{"location":"man/autodiff/","page":"AutoDiff","title":"AutoDiff","text":"or accessed in kernels written for KernelAbstractions.jl like for example the power flow equations (here in polar form):","category":"page"},{"location":"man/autodiff/","page":"AutoDiff","title":"AutoDiff","text":"@kernel function residual_kernel!(F, v_m, v_a,\n ybus_re_nzval, ybus_re_colptr, ybus_re_rowval,\n ybus_im_nzval, ybus_im_colptr, ybus_im_rowval,\n pinj, qinj, pv, pq, nbus)\n\n npv = size(pv, 1)\n npq = size(pq, 1)\n\n i = @index(Global, Linear)\n # REAL PV: 1:npv\n # REAL PQ: (npv+1:npv+npq)\n # IMAG PQ: (npv+npq+1:npv+2npq)\n fr = (i <= npv) ? pv[i] : pq[i - npv]\n F[i] -= pinj[fr]\n if i > npv\n F[i + npq] -= qinj[fr]\n end\n @inbounds for c in ybus_re_colptr[fr]:ybus_re_colptr[fr+1]-1\n to = ybus_re_rowval[c]\n aij = v_a[fr] - v_a[to]\n coef_cos = v_m[fr]*v_m[to]*ybus_re_nzval[c]\n coef_sin = v_m[fr]*v_m[to]*ybus_im_nzval[c]\n cos_val = cos(aij)\n sin_val = sin(aij)\n F[i] += coef_cos * cos_val + coef_sin * sin_val\n if i > npv\n F[npq + i] += coef_cos * sin_val - coef_sin * cos_val\n end\n end\nend","category":"page"},{"location":"man/autodiff/","page":"AutoDiff","title":"AutoDiff","text":"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","category":"page"},{"location":"man/autodiff/","page":"AutoDiff","title":"AutoDiff","text":"(Image: SIMD AD for point-wise vector product \\label{fig:simd})","category":"page"},{"location":"man/autodiff/","page":"AutoDiff","title":"AutoDiff","text":"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.","category":"page"},{"location":"man/autodiff/","page":"AutoDiff","title":"AutoDiff","text":"[1]: Griewank, Andreas, and Andrea Walther. Evaluating derivatives: principles and techniques of algorithmic differentiation. Society for Industrial and Applied Mathematics, 2008.","category":"page"}] +[{"location":"lib/formulations/","page":"Polar formulation","title":"Polar formulation","text":"CurrentModule = ExaPF\nconst PS = ExaPF.PowerSystem\nDocTestSetup = quote\n using ExaPF\nend","category":"page"},{"location":"lib/formulations/#Polar-formulation","page":"Polar formulation","title":"Polar formulation","text":"","category":"section"},{"location":"lib/formulations/#Generic-templates","page":"Polar formulation","title":"Generic templates","text":"","category":"section"},{"location":"lib/formulations/","page":"Polar formulation","title":"Polar formulation","text":"AbstractVariable\nAbstractFormulation\nState\nControl\n","category":"page"},{"location":"lib/formulations/#ExaPF.AbstractVariable","page":"Polar formulation","title":"ExaPF.AbstractVariable","text":"AbstractVariable\n\nVariables corresponding to a particular formulation.\n\n\n\n\n\n","category":"type"},{"location":"lib/formulations/#ExaPF.AbstractFormulation","page":"Polar formulation","title":"ExaPF.AbstractFormulation","text":"AbstractFormulation\n\nInterface between the data and the mathemical formulation.\n\n\n\n\n\n","category":"type"},{"location":"lib/formulations/#ExaPF.State","page":"Polar formulation","title":"ExaPF.State","text":"State <: AbstractVariable\n\nAll variables x depending on the variables Control u through the non-linear equation g(x u) = 0.\n\n\n\n\n\n","category":"type"},{"location":"lib/formulations/#ExaPF.Control","page":"Polar formulation","title":"ExaPF.Control","text":"Control <: AbstractVariable\n\nIndependent variables u used in the reduced-space formulation.\n\n\n\n\n\n","category":"type"},{"location":"lib/formulations/#Structure-and-variables","page":"Polar formulation","title":"Structure and variables","text":"","category":"section"},{"location":"lib/formulations/","page":"Polar formulation","title":"Polar formulation","text":"PolarForm\nBlockPolarForm\nload_polar\nNetworkStack\ninit!\n","category":"page"},{"location":"lib/formulations/#ExaPF.PolarForm","page":"Polar formulation","title":"ExaPF.PolarForm","text":"PolarForm{T, IT, VT, MT} <: AbstractPolarFormulation\n\nImplement the polar formulation associated to the network's equations.\n\nWrap a PS.PowerNetwork network to load the data on the target device (CPU() and CUDABackend() are currently supported).\n\nExample\n\njulia> const PS = ExaPF.PowerSystem;\n\njulia> network_data = PS.load_case(\"case9.m\");\n\njulia> polar = PolarForm(network_data, ExaPF.CPU())\nPolar formulation (instantiated on device CPU(false))\nNetwork characteristics:\n #buses: 9 (#slack: 1 #PV: 2 #PQ: 6)\n #generators: 3\n #lines: 9\ngiving a mathematical formulation with:\n #controls: 5\n #states : 14\n\n\n\n\n\n\n","category":"type"},{"location":"lib/formulations/#ExaPF.BlockPolarForm","page":"Polar formulation","title":"ExaPF.BlockPolarForm","text":"BlockPolarForm{T, IT, VT, MT} <: AbstractFormulation\n\nBlock polar formulation: duplicates k different polar models to evaluate them in parallel.\n\n\n\n\n\n","category":"type"},{"location":"lib/formulations/#ExaPF.load_polar","page":"Polar formulation","title":"ExaPF.load_polar","text":"load_polar(case, device=CPU(); dir=PS.EXADATA)\n\nLoad 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).\n\nExamples\n\njulia> polar = ExaPF.load_polar(\"case9\")\nPolar formulation (instantiated on device CPU(false))\nNetwork characteristics:\n #buses: 9 (#slack: 1 #PV: 2 #PQ: 6)\n #generators: 3\n #lines: 9\ngiving a mathematical formulation with:\n #controls: 5\n #states : 14\n\n\n\n\n\n\n","category":"function"},{"location":"lib/formulations/#ExaPF.NetworkStack","page":"Polar formulation","title":"ExaPF.NetworkStack","text":"NetworkStack{VT,VD,MT} <: AbstractNetworkStack{VT}\nNetworkStack(polar::PolarForm)\nNetworkStack(nbus::Int, ngen::Int, nlines::Int, VT::Type)\n\nStore the variables associated to the polar formulation. The variables are stored in the field input, ordered as follows\n\n input = [vmag ; vang ; pgen]\n\nThe object stores also intermediate variables needed in the expression tree, such as the LKMR basis ψ.\n\nNotes\n\nThe NetworkStack can be instantiated on the host or on the target device.\n\nExamples\n\njulia> polar = ExaPF.load_polar(\"case9\");\n\njulia> stack = ExaPF.NetworkStack(polar)\n21-elements NetworkStack{Vector{Float64}}\n\njulia> stack.vmag\n9-element Vector{Float64}:\n 1.0\n 1.0\n 1.0\n 1.0\n 1.0\n 1.0\n 1.0\n 1.0\n 1.0\n\n\n\n\n\n\n","category":"type"},{"location":"lib/formulations/#ExaPF.init!","page":"Polar formulation","title":"ExaPF.init!","text":"init!(polar::PolarForm, stack::NetworkStack)\n\nSet stack.input with the initial values specified in the base PS.PowerNetwork object.\n\n\n\n\n\n","category":"function"},{"location":"lib/formulations/","page":"Polar formulation","title":"Polar formulation","text":"The state and the control are defined as mapping:","category":"page"},{"location":"lib/formulations/","page":"Polar formulation","title":"Polar formulation","text":"mapping\n","category":"page"},{"location":"lib/formulations/#ExaPF.mapping","page":"Polar formulation","title":"ExaPF.mapping","text":"mapping(polar::PolarForm, ::Control)\n\nReturn the mapping associated to the Control() in NetworkStack according to the polar formulation PolarForm.\n\nExamples\n\njulia> polar = ExaPF.load_polar(\"case9\");\n\njulia> mapu = ExaPF.mapping(polar, Control())\n5-element Vector{Int64}:\n 1\n 2\n 3\n 20\n 21\n\n\n\n\n\n\nmapping(polar::PolarForm, ::State)\n\nReturn the mapping associated to the State() in NetworkStack according to the polar formulation PolarForm.\n\nExamples\n\njulia> polar = ExaPF.load_polar(\"case9\");\n\njulia> mapu = ExaPF.mapping(polar, State())\n14-element Vector{Int64}:\n 11\n 12\n 13\n 14\n 15\n 16\n 17\n 18\n 4\n 5\n 6\n 7\n 8\n 9\n\n\n\n\n\n\n","category":"function"},{"location":"lib/formulations/#Powerflow-solver","page":"Polar formulation","title":"Powerflow solver","text":"","category":"section"},{"location":"lib/formulations/","page":"Polar formulation","title":"Polar formulation","text":"run_pf\nnlsolve!\nNewtonRaphson\n","category":"page"},{"location":"lib/formulations/#ExaPF.run_pf","page":"Polar formulation","title":"ExaPF.run_pf","text":"run_pf(\n polar::PolarForm, stack::NetworkStack;\n rtol=1e-8, max_iter=20, verbose=0,\n)\n\nSolve 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.\n\nThe algorithm stops when a tolerance rtol or a maximum number of iterations maxiter is reached.\n\nArguments\n\npolar::AbstractFormulation: formulation of the power flow equation\nstack::NetworkStack: initial values in the network\n\nExamples\n\njulia> polar = ExaPF.load_polar(\"case9\");\n\njulia> stack = ExaPF.NetworkStack(polar);\n\njulia> conv = run_pf(polar, stack);\n\njulia> conv.has_converged\ntrue\n\njulia> conv.n_iterations\n4\n\n\n\n\n\n","category":"function"},{"location":"lib/formulations/#ExaPF.nlsolve!","page":"Polar formulation","title":"ExaPF.nlsolve!","text":"nlsolve!(\n algo::NewtonRaphson,\n jac::Jacobian,\n stack::NetworkStack;\n linear_solver=DirectSolver(jac.J),\n nl_buffer=NLBuffer(size(jac, 2)),\n)\n\nSolve the nonlinear system of equations g(x) = 0 with a NewtonRaphson algorithm. At each iteration, we update the variable x as\n\n x_k+1 = x_k - (g_k)^-1 g(x_k)\n\n\ntill g(x_k) ε_tol\n\nIn the implementation,\n\nthe function g is specified in jac.func,\nthe initial variable x_0 in stack::NetworkStack (with mapping jac.map),\nthe Jacobian g is computed automatically in jac, with automatic differentiation.\n\nNote that stack is modified inplace during the iterations of algorithm.\n\nThe Jacobian jac should be instantied before calling this function. By default, the linear system (g_k)^-1 g(x_k) is solved using a LU factorization. You can specify a different linear solver by changing the optional argument linear_solver.\n\nArguments\n\nalgo::NewtonRaphon: Newton-Raphson object, storing the options of the algorithm\njac::Jacobian: Stores the function g and its Jacobian g. The Jacobian is updated with automatic differentiation.\nstack::NetworkStack: initial values\nlinear_solver::AbstractLinearSolver: linear solver used to compute the Newton step\nnl_buffer::NLBuffer: buffer storing the residual vector and the descent direction Δx. Can be reused to avoid unecessary allocations.\n\nExamples\n\njulia> polar = ExaPF.load_polar(\"case9\");\n\njulia> powerflow = ExaPF.PowerFlowBalance(polar) ∘ ExaPF.PolarBasis(polar);\n\njulia> jx = ExaPF.Jacobian(polar, powerflow, State());\n\njulia> stack = ExaPF.NetworkStack(polar);\n\njulia> conv = ExaPF.nlsolve!(NewtonRaphson(), jx, stack);\n\njulia> conv.has_converged\ntrue\n\njulia> conv.n_iterations\n4\n\n\n\n\n\n","category":"function"},{"location":"lib/formulations/#ExaPF.NewtonRaphson","page":"Polar formulation","title":"ExaPF.NewtonRaphson","text":"NewtonRaphson <: AbstractNonLinearSolver\n\nNewton-Raphson algorithm.\n\nAttributes\n\nmaxiter::Int (default 20): maximum number of iterations\ntol::Float64 (default 1e-8): tolerance of the algorithm\nverbose::Int (default 0): verbosity level\n\n\n\n\n\n","category":"type"},{"location":"lib/formulations/#Constraints","page":"Polar formulation","title":"Constraints","text":"","category":"section"},{"location":"lib/formulations/","page":"Polar formulation","title":"Polar formulation","text":"The different parts of the polar formulation are implemented in the following AbstractExpression:","category":"page"},{"location":"lib/formulations/","page":"Polar formulation","title":"Polar formulation","text":"PolarBasis\nPowerFlowBalance\nVoltageMagnitudeBounds\nPowerGenerationBounds\nLineFlows\n","category":"page"},{"location":"lib/formulations/#ExaPF.PolarBasis","page":"Polar formulation","title":"ExaPF.PolarBasis","text":"PolarBasis{VI, MT} <: AbstractExpression\nPolarBasis(polar::AbstractPolarFormulation)\n\nImplement the LKMR nonlinear basis. Takes as input the voltage magnitudes vmag and the voltage angles vang and returns\n\n beginaligned\n psi_ell^C(v theta) = v^f v^t cos(theta_f - theta_t) quad forall ell = 1 cdots n_ell \n psi_ell^S(v theta) = v^f v^t sin(theta_f - theta_t) quad forall ell = 1 cdots n_ell \n psi_k(v theta) = v_k^2 quad forall k = 1 cdots n_b\n endaligned\n\nDimension: 2 * n_lines + n_bus\n\nComplexity\n\n3 n_lines + n_bus mul, n_lines cos and n_lines sin\n\nExamples\n\njulia> polar = ExaPF.load_polar(\"case9\");\n\njulia> stack = ExaPF.NetworkStack(polar);\n\njulia> basis = ExaPF.PolarBasis(polar)\nPolarBasis (AbstractExpression)\n\njulia> basis(stack)\n27-element Vector{Float64}:\n 1.0\n 1.0\n 1.0\n 1.0\n 1.0\n 1.0\n 1.0\n 1.0\n 1.0\n 0.0\n ⋮\n 1.0\n 1.0\n 1.0\n 1.0\n 1.0\n 1.0\n 1.0\n 1.0\n 1.0\n\n\n\n\n\n\n","category":"type"},{"location":"lib/formulations/#ExaPF.PowerFlowBalance","page":"Polar formulation","title":"ExaPF.PowerFlowBalance","text":"PowerFlowBalance{VT, MT}\nPowerFlowBalance(polar)\n\nImplement 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:\n\nbeginaligned\n p_i = v_i sum_j^n v_j (g_ijcos(theta_i - theta_j) + b_ijsin(theta_i - theta_j)) \n i PV PQ \n q_i = v_i sum_j^n v_j (g_ijsin(theta_i - theta_j) - b_ijcos(theta_i - theta_j)) \n i PQ\nendaligned\n\nRequire composition with PolarBasis.\n\nDimension: n_pv + 2 * n_pq\n\nComplexity\n\n2 SpMV\n\nExamples\n\njulia> polar = ExaPF.load_polar(\"case9\");\n\njulia> stack = ExaPF.NetworkStack(polar);\n\njulia> powerflow = ExaPF.PowerFlowBalance(polar) ∘ ExaPF.PolarBasis(polar);\n\njulia> round.(powerflow(stack); digits=6)\n14-element Vector{Float64}:\n -1.63\n -0.85\n 0.0\n 0.9\n 0.0\n 1.0\n 0.0\n 1.25\n -0.167\n 0.042\n -0.2835\n 0.171\n -0.2275\n 0.259\n\njulia> run_pf(polar, stack); # solve powerflow equations\n\njulia> isapprox(powerflow(stack), zeros(14); atol=1e-8)\ntrue\n\n\n\n\n\n\n","category":"type"},{"location":"lib/formulations/#ExaPF.VoltageMagnitudeBounds","page":"Polar formulation","title":"ExaPF.VoltageMagnitudeBounds","text":"VoltageMagnitudeBounds\n\nImplement 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:\n\nv_pq^ v_pq v_pq^ \n\nDimension: n_pq\n\nComplexity\n\n1 copyto\n\nNote\n\nIn the reduced space, the constraints on the voltage magnitudes at PV nodes v_pv are taken into account when bounding the control u.\n\nExamples\n\njulia> polar = ExaPF.load_polar(\"case9\");\n\njulia> stack = ExaPF.NetworkStack(polar);\n\njulia> voltage_pq = ExaPF.VoltageMagnitudeBounds(polar);\n\njulia> voltage_pq(stack)\n6-element Vector{Float64}:\n 1.0\n 1.0\n 1.0\n 1.0\n 1.0\n 1.0\n\n\n\n\n\n\n","category":"type"},{"location":"lib/formulations/#ExaPF.PowerGenerationBounds","page":"Polar formulation","title":"ExaPF.PowerGenerationBounds","text":"PowerGenerationBounds{VT, MT}\nPowerGenerationBounds(polar)\n\nConstraints 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\n\np_gref^ p_gref p_gref^ \nC_g q_g^ C_g q_g C_g q_g^ \n\nRequire composition with PolarBasis.\n\nDimension: n_pv + 2 n_ref\n\nComplexity\n\n1 copyto, 1 SpMV\n\nExamples\n\njulia> polar = ExaPF.load_polar(\"case9\");\n\njulia> stack = ExaPF.NetworkStack(polar);\n\njulia> run_pf(polar, stack); # solve powerflow equations\n\njulia> power_generators = ExaPF.PowerGenerationBounds(polar) ∘ ExaPF.PolarBasis(polar);\n\njulia> round.(power_generators(stack); digits=6)\n4-element Vector{Float64}:\n 0.719547\n 0.24069\n 0.144601\n -0.03649\n\n\n\n\n\n\n","category":"type"},{"location":"lib/formulations/#ExaPF.LineFlows","page":"Polar formulation","title":"ExaPF.LineFlows","text":"LineFlows{VT, MT}\nLineFlows(polar)\n\nImplement thermal limit constraints on the lines of the network.\n\nRequire composition with PolarBasis.\n\nDimension: 2 * n_lines\n\nComplexity\n\n4 SpMV, 4 * n_lines quadratic, 2 * n_lines add\n\nExamples\n\njulia> polar = ExaPF.load_polar(\"case9\");\n\njulia> stack = ExaPF.NetworkStack(polar);\n\njulia> run_pf(polar, stack); # solve powerflow equations\n\njulia> line_flows = ExaPF.LineFlows(polar) ∘ ExaPF.PolarBasis(polar);\n\njulia> round.(line_flows(stack); digits=6)\n18-element Vector{Float64}:\n 0.575679\n 0.094457\n 0.379983\n 0.723832\n 0.060169\n 0.588673\n 2.657418\n 0.748943\n 0.295351\n 0.560817\n 0.112095\n 0.38625\n 0.728726\n 0.117191\n 0.585164\n 2.67781\n 0.726668\n 0.215497\n\n\n\n\n\n\n","category":"type"},{"location":"lib/formulations/#Objective","page":"Polar formulation","title":"Objective","text":"","category":"section"},{"location":"lib/formulations/","page":"Polar formulation","title":"Polar formulation","text":"The production costs is given in the AbstractExpression CostFunction:","category":"page"},{"location":"lib/formulations/","page":"Polar formulation","title":"Polar formulation","text":"CostFunction","category":"page"},{"location":"lib/formulations/#ExaPF.CostFunction","page":"Polar formulation","title":"ExaPF.CostFunction","text":"CostFunction{VT, MT} <: AutoDiff.AbstractExpression\nCostFunction(polar)\n\nImplement the quadratic cost function for OPF\n\n _g=1^n_g c_2g p_g^2 + c_1g p_g + c_0g\n\nRequire composition with PolarBasis to evaluate the cost of the reference generator.\n\nDimension: 1\n\nComplexity\n\n1 SpMV, 1 sum\n\nExamples\n\njulia> polar = ExaPF.load_polar(\"case9\");\n\njulia> stack = ExaPF.NetworkStack(polar);\n\njulia> cost = ExaPF.CostFunction(polar) ∘ ExaPF.PolarBasis(polar);\n\njulia> cost(stack)\n1-element Vector{Float64}:\n 4509.0275\n\n\n\n\n\n\n","category":"type"},{"location":"lib/formulations/#Composition-of-expressions","page":"Polar formulation","title":"Composition of expressions","text":"","category":"section"},{"location":"lib/formulations/","page":"Polar formulation","title":"Polar formulation","text":"The different expressions can be combined together in several different ways.","category":"page"},{"location":"lib/formulations/","page":"Polar formulation","title":"Polar formulation","text":"MultiExpressions\nComposedExpressions","category":"page"},{"location":"lib/formulations/#ExaPF.MultiExpressions","page":"Polar formulation","title":"ExaPF.MultiExpressions","text":"MultiExpressions <: AbstractExpression\n\nImplement expressions concatenation. Takes as input a vector of expressions [expr1,...,exprN] and concatenate them in a single expression mexpr, such that\n\n mexpr(x) = [expr1(x) ; expr2(x) ; ... ; exprN(x)]\n\n\n\n\n\n\n","category":"type"},{"location":"lib/formulations/#ExaPF.ComposedExpressions","page":"Polar formulation","title":"ExaPF.ComposedExpressions","text":"ComposedExpressions{Expr1<:PolarBasis, Expr2} <: AbstractExpression\n\nImplement expression composition. Takes as input two expressions expr1 and expr2 and returns a composed expression cexpr such that ``` cexpr(x) = expr2 ∘ expr1(x)\n\nNotes\n\nCurrently, only PolarBasis is supported for expr1.\n\n\n\n\n\n","category":"type"},{"location":"lib/autodiff/","page":"AutoDiff","title":"AutoDiff","text":"CurrentModule = ExaPF.AutoDiff","category":"page"},{"location":"lib/autodiff/#AutoDiffAPI","page":"AutoDiff","title":"AutoDiff","text":"","category":"section"},{"location":"lib/autodiff/#Variables","page":"AutoDiff","title":"Variables","text":"","category":"section"},{"location":"lib/autodiff/","page":"AutoDiff","title":"AutoDiff","text":"AbstractStack\n","category":"page"},{"location":"lib/autodiff/#ExaPF.AutoDiff.AbstractStack","page":"AutoDiff","title":"ExaPF.AutoDiff.AbstractStack","text":"AbstractStack{VT}\n\nAbstract variable storing the inputs and the intermediate values in the expression tree.\n\n\n\n\n\n","category":"type"},{"location":"lib/autodiff/#Expressions","page":"AutoDiff","title":"Expressions","text":"","category":"section"},{"location":"lib/autodiff/","page":"AutoDiff","title":"AutoDiff","text":"AbstractExpression\nadjoint!\n","category":"page"},{"location":"lib/autodiff/#ExaPF.AutoDiff.AbstractExpression","page":"AutoDiff","title":"ExaPF.AutoDiff.AbstractExpression","text":"AbstractExpression\n\nAbstract 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).\n\nForward mode\n\nThe direct evaluation of the function f is implemented as\n\n(expr::AbstractExpression)(output::VT, stack::AbstractStack{VT}) where VT<:AbstractArray\n\n\nthe input being specified in stack, the results being stored in the array output.\n\nReverse mode\n\nThe adjoint of the function is specified by the function adjoint!, with the signature:\n\nadjoint!(expr::AbstractExpression, ∂stack::AbstractStack{VT}, stack::AbstractStack{VT}, ̄v::VT) where VT<:AbstractArray\n\n\nThe variable stack stores the result of the direct evaluation, and is not modified in adjoint!. The results are stored inside the adjoint stack ∂stack.\n\n\n\n\n\n","category":"type"},{"location":"lib/autodiff/#ExaPF.AutoDiff.adjoint!","page":"AutoDiff","title":"ExaPF.AutoDiff.adjoint!","text":"adjoint!(expr::AbstractExpression, ∂stack::AbstractStack{VT}, stack::AbstractStack{VT}, ̄v::VT) where VT<:AbstractArray\n\nCompute 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!.\n\n\n\n\n\n","category":"function"},{"location":"lib/autodiff/#First-and-second-order-derivatives","page":"AutoDiff","title":"First and second-order derivatives","text":"","category":"section"},{"location":"lib/autodiff/","page":"AutoDiff","title":"AutoDiff","text":"AbstractJacobian\nAbstractHessianProd\nAbstractFullHessian","category":"page"},{"location":"lib/autodiff/#ExaPF.AutoDiff.AbstractJacobian","page":"AutoDiff","title":"ExaPF.AutoDiff.AbstractJacobian","text":"AbstractJacobian\n\nAutomatic differentiation for the compressed Jacobian of any nonlinear constraint h(x).\n\n\n\n\n\n","category":"type"},{"location":"lib/autodiff/#ExaPF.AutoDiff.AbstractHessianProd","page":"AutoDiff","title":"ExaPF.AutoDiff.AbstractHessianProd","text":"AbstractHessianProd\n\nReturns the adjoint-Hessian-vector product λ^ H v of any nonlinear constraint h(x).\n\n\n\n\n\n","category":"type"},{"location":"lib/autodiff/#ExaPF.AutoDiff.AbstractFullHessian","page":"AutoDiff","title":"ExaPF.AutoDiff.AbstractFullHessian","text":"AbstractHessianProd\n\nFull sparse Hessian H of any nonlinear constraint h(x).\n\n\n\n\n\n","category":"type"},{"location":"lib/autodiff/#Utils","page":"AutoDiff","title":"Utils","text":"","category":"section"},{"location":"lib/autodiff/","page":"AutoDiff","title":"AutoDiff","text":"seed!\nseed_coloring!\npartials!\nset_value!","category":"page"},{"location":"lib/autodiff/#ExaPF.AutoDiff.seed!","page":"AutoDiff","title":"ExaPF.AutoDiff.seed!","text":"seed!(\n H::AbstractHessianProd,\n v::AbstractVector{T},\n) where {T}\n\nSeed the duals with v to compute the Hessian vector product λ^ H v.\n\n\n\n\n\n","category":"function"},{"location":"lib/autodiff/#ExaPF.AutoDiff.seed_coloring!","page":"AutoDiff","title":"ExaPF.AutoDiff.seed_coloring!","text":"seed_coloring!(\n M::Union{AbstractJacobian, AbstractFullHessian}\n coloring::AbstractVector,\n)\n\nSeed the duals with the coloring based seeds to compute the Jacobian or Hessian M.\n\n\n\n\n\n","category":"function"},{"location":"lib/autodiff/#ExaPF.AutoDiff.partials!","page":"AutoDiff","title":"ExaPF.AutoDiff.partials!","text":"partials!(jac::AbstractJacobian)\n\nExtract partials from Jacobian jac in jac.J.\n\n\n\n\n\npartials!(hess::AbstractFullHessian)\n\nExtract partials from Hessian hess into hess.H.\n\n\n\n\n\n","category":"function"},{"location":"lib/autodiff/#ExaPF.AutoDiff.set_value!","page":"AutoDiff","title":"ExaPF.AutoDiff.set_value!","text":"set_value!(\n jac,\n primals::AbstractVector{T}\n) where {T}\n\nSet values of ForwardDiff.Dual numbers in jac to primals.\n\n\n\n\n\n","category":"function"},{"location":"man/powersystem/","page":"PowerSystem","title":"PowerSystem","text":"CurrentModule = ExaPF\nDocTestSetup = quote\n using ExaPF\n const PS = ExaPF.PowerSystem\nend\nDocTestFilters = [r\"ExaPF\"]","category":"page"},{"location":"man/powersystem/#PowerSystem","page":"PowerSystem","title":"PowerSystem","text":"","category":"section"},{"location":"man/powersystem/","page":"PowerSystem","title":"PowerSystem","text":"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.","category":"page"},{"location":"man/powersystem/","page":"PowerSystem","title":"PowerSystem","text":"We usually load the PowerSystem system submodule with the alias PS:","category":"page"},{"location":"man/powersystem/","page":"PowerSystem","title":"PowerSystem","text":"julia> PS = ExaPF.PowerSystem\n","category":"page"},{"location":"man/powersystem/#Description","page":"PowerSystem","title":"Description","text":"","category":"section"},{"location":"man/powersystem/","page":"PowerSystem","title":"PowerSystem","text":"The electrical power system is represented as a linear, lumped network which has to satisfy the Kirchhoff laws:","category":"page"},{"location":"man/powersystem/","page":"PowerSystem","title":"PowerSystem","text":" bmi = bmYbmv ","category":"page"},{"location":"man/powersystem/","page":"PowerSystem","title":"PowerSystem","text":"where bmi bmv in mathbbC^N_B are the current and voltage vectors associated to the system and bmY in mathbbC^N_B times N_B is the admittance matrix. These equations are often rewritten in terms of apparent powers:","category":"page"},{"location":"man/powersystem/","page":"PowerSystem","title":"PowerSystem","text":" bms = bmp + jbmq = textitdiag(bmv^*) bmYbmv","category":"page"},{"location":"man/powersystem/","page":"PowerSystem","title":"PowerSystem","text":"Using polar representation of the voltage vector, such as bmv = ve^j theta, each bus i=1 cdots N_B must satisfy the power balance equations:","category":"page"},{"location":"man/powersystem/","page":"PowerSystem","title":"PowerSystem","text":"beginaligned\n p_i = v_i sum_j^n v_j (g_ijcos(theta_i - theta_j) + b_ijsin(theta_i - theta_j)) \n q_i = v_i sum_j^n v_j (g_ijsin(theta_i - theta_j) - b_ijcos(theta_i - theta_j)) \nendaligned","category":"page"},{"location":"man/powersystem/","page":"PowerSystem","title":"PowerSystem","text":"where each bus i has variables p_i q_i v_i theta_i and the topology of the network is defined by a non-negative value of the admittance between two buses i and j, y_ij = g_ij + ib_ij.","category":"page"},{"location":"man/powersystem/#The-PowerNetwork-Object","page":"PowerSystem","title":"The PowerNetwork Object","text":"","category":"section"},{"location":"man/powersystem/","page":"PowerSystem","title":"PowerSystem","text":"Currently we can create a PS.PowerNetwork object by parsing a MATPOWER data file.","category":"page"},{"location":"man/powersystem/","page":"PowerSystem","title":"PowerSystem","text":"julia> datafile = \"case9.m\";\n\njulia> ps = PS.load_case(datafile)\nPowerNetwork object with:\n Buses: 9 (Slack: 1. PV: 2. PQ: 6)\n Generators: 3.\n","category":"page"},{"location":"man/powersystem/","page":"PowerSystem","title":"PowerSystem","text":"Then, using multiple dispatch, we have defined a set of abstract data types and getter functions which allow us to retrieve information from the PowerNetwork object","category":"page"},{"location":"man/powersystem/","page":"PowerSystem","title":"PowerSystem","text":"julia> PS.get(ps, PS.NumberOfPQBuses())\n6\n\njulia> PS.get(ps, PS.NumberOfPVBuses())\n2\n\njulia> PS.get(ps, PS.NumberOfSlackBuses())\n1","category":"page"},{"location":"man/linearsolver/","page":"Linear Solvers","title":"Linear Solvers","text":"CurrentModule = ExaPF\nDocTestSetup = quote\n using ExaPF\n const Precondition = ExaPF.Precondition\n const Iterative = ExaPF.Iterative\nend\nDocTestFilters = [r\"ExaPF\"]","category":"page"},{"location":"man/linearsolver/#Linear-Solver","page":"Linear Solvers","title":"Linear Solver","text":"","category":"section"},{"location":"man/linearsolver/#Overview","page":"Linear Solvers","title":"Overview","text":"","category":"section"},{"location":"man/linearsolver/","page":"Linear Solvers","title":"Linear Solvers","text":"As mentioned before, a linear solver is required to compute the Newton step in","category":"page"},{"location":"man/linearsolver/","page":"Linear Solvers","title":"Linear Solvers","text":"dx .= jacobian(x)\\f(x)","category":"page"},{"location":"man/linearsolver/","page":"Linear Solvers","title":"Linear Solvers","text":"Our package supports the following linear solvers:","category":"page"},{"location":"man/linearsolver/","page":"Linear Solvers","title":"Linear Solvers","text":"cuSOLVER with csrlsvqr (GPU),\nKrylov.jl with dqgmres and bicgstab (CPU/GPU),\nUMFPACK through the default Julia \\ operator (CPU),\ngeneric BiCGSTAB implementation [Vorst1992] (CPU/GPU),\nor any linear solver wrapped in LinearAlgebra.","category":"page"},{"location":"man/linearsolver/#Preconditioning","page":"Linear Solvers","title":"Preconditioning","text":"","category":"section"},{"location":"man/linearsolver/","page":"Linear Solvers","title":"Linear Solvers","text":"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.","category":"page"},{"location":"man/linearsolver/","page":"Linear Solvers","title":"Linear Solvers","text":"The block-Jacobi preconditioner used in ExaPF has been added to KrylovPreconditioners.jl","category":"page"},{"location":"man/linearsolver/","page":"Linear Solvers","title":"Linear Solvers","text":"The Jacobian is partitioned into a dense block diagonal structure using Metis.jl, where each block is inverted to build our preconditioner P.","category":"page"},{"location":"man/linearsolver/","page":"Linear Solvers","title":"Linear Solvers","text":"(Image: Dense block Jacobi preconditioner \\label{fig:preconditioner})","category":"page"},{"location":"man/linearsolver/","page":"Linear Solvers","title":"Linear Solvers","text":"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.","category":"page"},{"location":"man/linearsolver/","page":"Linear Solvers","title":"Linear Solvers","text":"CUDA.@sync pivot, info = CUDA.CUBLAS.getrf_batched!(blocks, true)\nCUDA.@sync pivot, info, p.cuJs = CUDA.CUBLAS.getri_batched(blocks, pivot)","category":"page"},{"location":"man/linearsolver/","page":"Linear Solvers","title":"Linear Solvers","text":"Assuming that other vendors will provide such batched BLAS APIs, this code is portable to other GPU architectures.","category":"page"},{"location":"man/linearsolver/","page":"Linear Solvers","title":"Linear Solvers","text":"[Vorst1992]: Vorst, H. A. van der. 1992. “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 (2): 631–44","category":"page"},{"location":"artifact/#ExaData-Artifact","page":"ExaData Artifact","title":"ExaData Artifact","text":"","category":"section"},{"location":"artifact/","page":"ExaData Artifact","title":"ExaData Artifact","text":"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.","category":"page"},{"location":"tutorials/direct_solver/","page":"Power flow: direct solver","title":"Power flow: direct solver","text":"CurrentModule = ExaPF\nDocTestSetup = quote\n using ExaPF\n const LS = ExaPF.LinearSolvers\nend\nDocTestFilters = [r\"ExaPF\"]","category":"page"},{"location":"tutorials/direct_solver/","page":"Power flow: direct solver","title":"Power flow: direct solver","text":"using ExaPF\nusing KLU\nusing LinearAlgebra\nconst LS = ExaPF.LinearSolvers\n","category":"page"},{"location":"tutorials/direct_solver/#Direct-solvers-for-power-flow","page":"Power flow: direct solver","title":"Direct solvers for power flow","text":"","category":"section"},{"location":"tutorials/direct_solver/","page":"Power flow: direct solver","title":"Power flow: direct solver","text":"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","category":"page"},{"location":"tutorials/direct_solver/","page":"Power flow: direct solver","title":"Power flow: direct solver","text":"g(x p) = 0","category":"page"},{"location":"tutorials/direct_solver/","page":"Power flow: direct solver","title":"Power flow: direct solver","text":"where g mathbbR^n_x times mathbbR^n_p to mathbbR^n_x is a nonlinear function encoding the power flow equations.","category":"page"},{"location":"tutorials/direct_solver/","page":"Power flow: direct solver","title":"Power flow: direct solver","text":"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","category":"page"},{"location":"tutorials/direct_solver/","page":"Power flow: direct solver","title":"Power flow: direct solver","text":"(nabla_x g_k) Delta x = - g_k","category":"page"},{"location":"tutorials/direct_solver/","page":"Power flow: direct solver","title":"Power flow: direct solver","text":"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.","category":"page"},{"location":"tutorials/direct_solver/","page":"Power flow: direct solver","title":"Power flow: direct solver","text":"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","category":"page"},{"location":"tutorials/direct_solver/","page":"Power flow: direct solver","title":"Power flow: direct solver","text":"P (nabla_x g) Q = LU","category":"page"},{"location":"tutorials/direct_solver/","page":"Power flow: direct solver","title":"Power flow: direct solver","text":"Once the matrix factorized, solving the linear system just translates to perform two backsolves with the triangular matrices L and U.","category":"page"},{"location":"tutorials/direct_solver/","page":"Power flow: direct solver","title":"Power flow: direct solver","text":"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.","category":"page"},{"location":"tutorials/direct_solver/#UMFPACK-(CPU,-default)","page":"Power flow: direct solver","title":"UMFPACK (CPU, default)","text":"","category":"section"},{"location":"tutorials/direct_solver/","page":"Power flow: direct solver","title":"Power flow: direct solver","text":"By default, ExaPF employs the linear solver UMFPACK to solve the linear system, as UMFPACK is shipped automatically in Julia.","category":"page"},{"location":"tutorials/direct_solver/","page":"Power flow: direct solver","title":"Power flow: direct solver","text":"In the LinearSolvers submodule, this is how the wrapper is implemented:","category":"page"},{"location":"tutorials/direct_solver/","page":"Power flow: direct solver","title":"Power flow: direct solver","text":"struct DirectSolver{Fac} <: AbstractLinearSolver\n factorization::Fac\nend\nDirectSolver(J::AbstractMatrix) = DirectSolver(lu(J))\n","category":"page"},{"location":"tutorials/direct_solver/","page":"Power flow: direct solver","title":"Power flow: direct solver","text":"By default, the constructor takes as input the initial Jacobian J and factorizes it by calling lu(J), which in Julia translates to a factorization with UMFPACK. Then, inside the function nlsolve! we refactorize the matrix at each iteration by calling the function LinearSolvers.update!","category":"page"},{"location":"tutorials/direct_solver/","page":"Power flow: direct solver","title":"Power flow: direct solver","text":"function update!(s::DirectSolver, J::AbstractMatrix)\n LinearAlgebra.lu!(s.factorization, J) # Update factorization inplace\nend","category":"page"},{"location":"tutorials/direct_solver/","page":"Power flow: direct solver","title":"Power flow: direct solver","text":"This function uses the function LinearAlgebra.lu! to update the factorization inplace. The backsolve is computed by calling the LinearAlgebra.ldiv! function:","category":"page"},{"location":"tutorials/direct_solver/","page":"Power flow: direct solver","title":"Power flow: direct solver","text":"function ldiv!(s::DirectSolver, y::AbstractVector, J::AbstractMatrix, x::AbstractVector)\n LinearAlgebra.ldiv!(y, s.factorization, x) # Forward-backward solve\n return 0\nend","category":"page"},{"location":"tutorials/direct_solver/","page":"Power flow: direct solver","title":"Power flow: direct solver","text":"We notice that the code has been designed to support any factorization routines implementing the two routines LinearAlgebra.lu! and LinearAlgebra.ldiv!.","category":"page"},{"location":"tutorials/direct_solver/","page":"Power flow: direct solver","title":"Power flow: direct solver","text":"Before comparing with other linear solvers, we solve a large scale power flow instance with UMFPACK to give us a reference.","category":"page"},{"location":"tutorials/direct_solver/","page":"Power flow: direct solver","title":"Power flow: direct solver","text":"polar = ExaPF.load_polar(\"case9241pegase.m\")\nstack = ExaPF.NetworkStack(polar)\npf_solver = NewtonRaphson(tol=1e-10, verbose=2) # power flow solver\nfunc = ExaPF.PowerFlowBalance(polar) ∘ ExaPF.PolarBasis(polar) # power flow func\njx = ExaPF.Jacobian(polar, func, State()) # init AD\nExaPF.nlsolve!(pf_solver, jx, stack)","category":"page"},{"location":"tutorials/direct_solver/#KLU-(CPU)","page":"Power flow: direct solver","title":"KLU (CPU)","text":"","category":"section"},{"location":"tutorials/direct_solver/","page":"Power flow: direct solver","title":"Power flow: direct solver","text":"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!","category":"page"},{"location":"tutorials/direct_solver/","page":"Power flow: direct solver","title":"Power flow: direct solver","text":"LinearAlgebra.lu!(K::KLU.KLUFactorization, J) = KLU.klu!(K, J)\n","category":"page"},{"location":"tutorials/direct_solver/","page":"Power flow: direct solver","title":"Power flow: direct solver","text":"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:","category":"page"},{"location":"tutorials/direct_solver/","page":"Power flow: direct solver","title":"Power flow: direct solver","text":"klu_factorization = KLU.klu(jx.J)\nklu_solver = LS.DirectSolver(klu_factorization)\n","category":"page"},{"location":"tutorials/direct_solver/","page":"Power flow: direct solver","title":"Power flow: direct solver","text":"and pass it to nlsolve!:","category":"page"},{"location":"tutorials/direct_solver/","page":"Power flow: direct solver","title":"Power flow: direct solver","text":"ExaPF.init!(polar, stack) # reinit stack\nExaPF.nlsolve!(pf_solver, jx, stack; linear_solver=klu_solver)\n","category":"page"},{"location":"tutorials/direct_solver/","page":"Power flow: direct solver","title":"Power flow: direct solver","text":"We observe KLU reduces considerably the time spent in the linear solver.","category":"page"},{"location":"tutorials/direct_solver/#cusolverRF-(CUDA)","page":"Power flow: direct solver","title":"cusolverRF (CUDA)","text":"","category":"section"},{"location":"tutorials/direct_solver/","page":"Power flow: direct solver","title":"Power flow: direct solver","text":"cusolverRF is an efficient LU refactorization routine implemented in CUDA. It is wrapped in Julia inside the package CUSOLVERRF.jl:","category":"page"},{"location":"tutorials/direct_solver/","page":"Power flow: direct solver","title":"Power flow: direct solver","text":"using CUSOLVERRF","category":"page"},{"location":"tutorials/direct_solver/","page":"Power flow: direct solver","title":"Power flow: direct solver","text":"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.","category":"page"},{"location":"tutorials/direct_solver/","page":"Power flow: direct solver","title":"Power flow: direct solver","text":"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.","category":"page"},{"location":"tutorials/direct_solver/","page":"Power flow: direct solver","title":"Power flow: direct solver","text":"CUSOLVERRF.jl follows the LinearAlgebra's interface, so we can use it directly in ExaPF. We first have to instantiate everything on the GPU:","category":"page"},{"location":"tutorials/direct_solver/","page":"Power flow: direct solver","title":"Power flow: direct solver","text":"using CUDA\npolar_gpu = ExaPF.load_polar(\"case9241pegase.m\", CUDABackend())\nstack_gpu = ExaPF.NetworkStack(polar_gpu)\nfunc_gpu = ExaPF.PowerFlowBalance(polar_gpu) ∘ ExaPF.PolarBasis(polar_gpu)\njx_gpu = ExaPF.Jacobian(polar_gpu, func_gpu, State()) # init AD","category":"page"},{"location":"tutorials/direct_solver/","page":"Power flow: direct solver","title":"Power flow: direct solver","text":"We can instantiate a new cusolverRF's instance as","category":"page"},{"location":"tutorials/direct_solver/","page":"Power flow: direct solver","title":"Power flow: direct solver","text":"rf_fac = CUSOLVERRF.RFLU(jx_gpu.J)\nrf_solver = LS.DirectSolver(rf_fac)\n","category":"page"},{"location":"tutorials/direct_solver/","page":"Power flow: direct solver","title":"Power flow: direct solver","text":"Then, we are able to solve the power flow entirely on the GPU, simply as","category":"page"},{"location":"tutorials/direct_solver/","page":"Power flow: direct solver","title":"Power flow: direct solver","text":"ExaPF.nlsolve!(pf_solver, jx_gpu, stack_gpu; linear_solver=rf_solver)\n","category":"page"},{"location":"man/benchmark/","page":"Benchmark","title":"Benchmark","text":"CurrentModule = ExaPF\nDocTestSetup = quote\n using ExaPF\nend\nDocTestFilters = [r\"ExaPF\"]","category":"page"},{"location":"man/benchmark/#Benchmark","page":"Benchmark","title":"Benchmark","text":"","category":"section"},{"location":"man/benchmark/","page":"Benchmark","title":"Benchmark","text":"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. Bicgstab), 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.","category":"page"},{"location":"man/benchmark/","page":"Benchmark","title":"Benchmark","text":"Running","category":"page"},{"location":"man/benchmark/","page":"Benchmark","title":"Benchmark","text":"julia --project benchmark/benchmarks.jl Bicgstab CUDABackend case300.m","category":"page"},{"location":"man/benchmark/","page":"Benchmark","title":"Benchmark","text":"yields","category":"page"},{"location":"man/benchmark/","page":"Benchmark","title":"Benchmark","text":"Bicgstab, CUDABackend, case300.m, 69.0, 3.57, 43.7, true","category":"page"},{"location":"man/benchmark/","page":"Benchmark","title":"Benchmark","text":"The first three fields are the settings of the benchmark run. They are followed by three timings in milliseconds:","category":"page"},{"location":"man/benchmark/","page":"Benchmark","title":"Benchmark","text":"the time taken by the Newton-Raphson algorithm to solve the power flow,\nthe timings for the Jacobian accumulation using AutoDiff,\nand the time for the linear solver, including the preconditioner.","category":"page"},{"location":"man/benchmark/","page":"Benchmark","title":"Benchmark","text":"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","category":"page"},{"location":"man/benchmark/","page":"Benchmark","title":"Benchmark","text":"powerflow_solver = NewtonRaphson(tol=ntol)","category":"page"},{"location":"man/benchmark/","page":"Benchmark","title":"Benchmark","text":"to one of the following options:","category":"page"},{"location":"man/benchmark/","page":"Benchmark","title":"Benchmark","text":"powerflow_solver = NewtonRaphson(tol=ntol, verbose=VERBOSE_LEVEL_NONE)\npowerflow_solver = NewtonRaphson(tol=ntol, verbose=VERBOSE_LEVEL_LOW)\npowerflow_solver = NewtonRaphson(tol=ntol, verbose=VERBOSE_LEVEL_MEDIUM)\npowerflow_solver = NewtonRaphson(tol=ntol, verbose=VERBOSE_LEVEL_HIGH)","category":"page"},{"location":"man/benchmark/","page":"Benchmark","title":"Benchmark","text":"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.","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"CurrentModule = ExaPF\nDocTestSetup = quote\n using ExaPF\n using LazyArtifacts\n import ExaPF: AutoDiff\nend","category":"page"},{"location":"man/formulations/#Formulations","page":"Polar formulation","title":"Formulations","text":"","category":"section"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"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.","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"In what follows, we denote by v in mathbbR^n_b the voltage magnitudes, theta in mathbbR^n_b the voltage angles and p_g q_g in mathbbR^n_g the active and reactive power generations. The active and reactive loads are denoted respectively by p_d q_d in mathbbR^n_b.","category":"page"},{"location":"man/formulations/#Power-flow-model","page":"Polar formulation","title":"Power flow model","text":"","category":"section"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"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 mathbbR^n_b times mathbbR^n_b to mathbbR^2n_ell + n_b. If we introduce the intermediate expressions","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":" psi_ell^C(v theta) = v^f v^t cos(theta_f - theta_t) quad forall ell = 1 cdots n_ell \n psi_ell^S(v theta) = v^f v^t sin(theta_f - theta_t) quad forall ell = 1 cdots n_ell \n psi_k(v theta) = v_k^2 quad forall k = 1 cdots n_b","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"the basis psi is defined as","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":" psi(v theta) = psi_ell^C(v theta)^top psi_ell^S(v theta)^top psi_k(v theta)^top ","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"The basis psi encodes all the nonlinearities in the problem. For instance, the power flow equations rewrite directly as","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":" beginbmatrix\n C_g p_g - p_d \n C_g q_g - q_d\n endbmatrix\n +\n underbrace\n beginbmatrix\n - hatG^c - hatB^s -G^d \n hatB^c - hatG^s B^d\n endbmatrix\n _M\n psi(v theta)\n = 0","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"with C_g in mathbbR^n_b times n_g the bus-generators incidence matrix, and the matrices B G defined from the admittance matrix Y_b of the network.","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"Similarly, the line flows rewrite","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":" beginbmatrix\n s_p^f s_q^f\n endbmatrix\n =\n overbrace\n beginbmatrix\n G_ft B_ft G_ff C_f^top \n -B_ft G_ft -B_ff C_f^top\n endbmatrix\n ^L_line^f\n psi(v theta) \n beginbmatrix\n s_p^t s_q^t\n endbmatrix\n =\n underbrace\n beginbmatrix\n G_tf B_tf G_tt C_t^top \n -B_tf G_tf -B_tt C_t^top\n endbmatrix\n _L_line^t\n psi(v theta)","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"with C_f in mathbbR^n_b times n_ell the bus-from incidence matrix and C_t in mathbbR^n_b times n_ell the bus-to incidence matrix. Then, the line flows constraints write directly with the quadratic expressions:","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":" (s_p^f)^2 + (s_q^f)^2 leq (s^max)^2 quad \n (s_p^t)^2 + (s_q^t)^2 leq (s^max)^2 quad ","category":"page"},{"location":"man/formulations/#Why-is-this-model-advantageous?","page":"Polar formulation","title":"Why is this model advantageous?","text":"","category":"section"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"Implementing the model LTMR2020 is not difficult once the basis function psi has been defined. Indeed, if we select a subset of the power flow equations (as usual, associated to the active injections at PV nodes, and active and reactive injections at PQ nodes), we get","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":" C_eq p_g + M_eq psi + tau = 0","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"with C_eq defined from the bus-generator incidence matrix C_g, M_eq a subset of the matrix M, tau a constant depending on the loads in the problem. Note that C_eq and M_eq are sparse matrices, so the expression can be implemented efficiently with sparse linear algebra operations (2 SpMV operations, 2 vector additions). The same holds true for the line flow constraints, evaluated with 2 SpMV operations:","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":" s^f = L_line^f psi quad\n s^t = L_line^t psi ","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"In ExaPF, all nonlinear expressions are written as linear operations depending on the nonlinear basis psi. By doing so, all the unstructured sparsity of the power flow problem is directly handled inside the sparse linear algebra library (cusparse on CUDA GPU, SuiteSparse on the CPU).","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"In what follows, we detail the implementation of the LTMR2020 model in ExaPF.","category":"page"},{"location":"man/formulations/#How-to-instantiate-the-inputs?","page":"Polar formulation","title":"How to instantiate the inputs?","text":"","category":"section"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"We have implemented the LTMR2020 model in ExaPF, both on the CPU and on CUDA GPU. All the operations have been rewritten in a vectorized fashion. Every model depends on inputs we propagate forward with functions. In ExaPF, the inputs will be specified in a NetworkStack <: AbstractStack. The functions will be implemented as AutoDiff.AbstractExpression.","category":"page"},{"location":"man/formulations/#Specifying-inputs-in-NetworkStack","page":"Polar formulation","title":"Specifying inputs in NetworkStack","text":"","category":"section"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"Our three inputs are (v theta p_g) in mathbbR^2n_b + n_g (voltage magnitude, voltage angle, power generations). The basis psi is considered as an intermediate expression.","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"We store all inputs in a NetworkStack structure:","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"struct NetworkStack{VT} <: AbstractStack\n input::VT\n vmag::VT # voltage magnitudes (view)\n vang::VT # voltage angles (view)\n pgen::VT # active power generations (view)\n ψ::VT # nonlinear basis ψ(vmag, vang)\nend","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"All the inputs are specified in the vector input. The three vectors vmag, vang and pgen are views porting on input, and are defined mostly for convenience. By convention the vector input is ordered as [vmag; vang; pgen]:","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"# Define dimension of the problem\njulia> nbus, ngen, nlines = 3, 2, 4;\n\njulia> stack = ExaPF.NetworkStack(nbus, ngen, nlines, 0, Vector{Float64}, Vector{Float64})\n8-elements NetworkStack{Vector{Float64}}\n\njulia> stack.input\n8-element Vector{Float64}:\n 0.0\n 0.0\n 0.0\n 0.0\n 0.0\n 0.0\n 0.0\n 0.0\n\njulia> stack.vmag .= 1;\n\njulia> stack.vang .= 2;\n\njulia> stack.pgen .= 3;\n\njulia> stack.input\n8-element Vector{Float64}:\n 1.0\n 1.0\n 1.0\n 2.0\n 2.0\n 2.0\n 3.0\n 3.0","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"The basis vector ψ is an intermediate expression, whose values depend on the inputs.","category":"page"},{"location":"man/formulations/#Defining-a-state-and-a-control","page":"Polar formulation","title":"Defining a state and a control","text":"","category":"section"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"In the reduced space method, we have to split the variables in a state x and a control u. By default, we define","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":" x = (theta_pv theta_pq v_pq) quad\n x = (v_ref v_pv p_ggenpv) ","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"and the control, and was not flexible. In the new implementation, we define the state and the control as two mappings porting on the vector stack.input (which itself stores all the inputs in the problem):","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"julia> nbus, ngen, nlines = 4, 3, 4;\n\njulia> stack = ExaPF.NetworkStack(nbus, ngen, nlines, 0, Vector{Float64}, Vector{Float64});\n\njulia> stack.input .= 1:length(stack.input); # index array input\n\njulia> ref, pv, pq, genpv = [1], [2], [3, 4], [2, 3];\n\njulia> mapx = [nbus .+ pv; nbus .+ pq; pq];\n\njulia> mapu = [ref; pv; 2*nbus .+ genpv];\n\njulia> x = @view stack.input[mapx]\n5-element view(::Vector{Float64}, [6, 7, 8, 3, 4]) with eltype Float64:\n 6.0\n 7.0\n 8.0\n 3.0\n 4.0\n\njulia> u = @view stack.input[mapu]\n4-element view(::Vector{Float64}, [1, 2, 10, 11]) with eltype Float64:\n 1.0\n 2.0\n 10.0\n 11.0","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"By doing so, the values of the state and the control are directly stored inside the NetworkStack structure, avoiding to duplicate values in the memory.","category":"page"},{"location":"man/formulations/#How-to-manipulate-the-expressions?","page":"Polar formulation","title":"How to manipulate the expressions?","text":"","category":"section"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"ExaPF implements the different functions required to implement the optimal power flow problem with the polar formulation:","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"PowerFlowBalance: power flow balance equations\nPowerGenerationBounds: bounds on the power generation\nLineFlows: line flow constraints\nCostFunction: quadratic cost function","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"Each function follows the LTMR2020 model and depends on the basis function psi(v theta), here implemented in the PolarBasis function.","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"We demonstrate how to use the different functions on the case9 instance. The procedure remains the same for all power network.","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"julia> polar = ExaPF.load_polar(\"case9.m\");\n\njulia> stack = ExaPF.NetworkStack(polar);\n","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"note: Note\nAll the code presented below is agnostic with regards to the specific device (CPU, CUDABackend...) we are using. By default, ExaPF computes the expressions on the CPU. Deporting the computation on a CUDABackend simply translates to instantiate the PolarForm structure on the GPU: polar = PolarForm(\"case9.m\", CUDABackend()).","category":"page"},{"location":"man/formulations/#Interface","page":"Polar formulation","title":"Interface","text":"","category":"section"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"All functions are following AutoDiff.AbstractExpression's interface. The structure of the network is specified by the PolarForm we pass as an argument in the constructor. For instance, we build a new PolarBasis expression associated to case9 directly as","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"julia> basis = ExaPF.PolarBasis(polar)\nPolarBasis (AbstractExpression)\n","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"Each expression as a given dimension, given by","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"julia> length(basis)\n27\n","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"In ExaPF, the inputs and the parameters are stored inside a NetworkStack structure. Evaluating the basis psi naturally translates to","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"julia> basis(stack)\n27-element Vector{Float64}:\n 1.0\n 1.0\n 1.0\n 1.0\n 1.0\n 1.0\n 1.0\n 1.0\n 1.0\n 0.0\n ⋮\n 1.0\n 1.0\n 1.0\n 1.0\n 1.0\n 1.0\n 1.0\n 1.0\n 1.0","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"This function call allocates a vector psi with 27 elements and evaluates the basis associated to the LTMR2020 model. To avoid unnecessary allocation, one can preallocate the vector psi:","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"julia> psi = zeros(length(basis)) ;\n\njulia> basis(psi, stack);\n","category":"page"},{"location":"man/formulations/#Compose-expressions-together","page":"Polar formulation","title":"Compose expressions together","text":"","category":"section"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"In the LTMR2020 model, the polar basis psi(v theta) depends only on the voltage magnitudes and the voltage angles. However, this is not the case for the other functions (power flow balance, line flows, ...), which all depends on the basis psi(v theta).","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"In ExaPF, one has to build manually the vectorized expression tree associated to the power flow model. Luckily, evaluating the LTMR2020 simply amounts to compose functions together with the polar basis psi(v theta). ExaPF overloads the function ∘ to compose functions with a PolarBasis instance. The power flow balance can be evaluated as","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"julia> pflow = ExaPF.PowerFlowBalance(polar) ∘ basis;\n","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"which returns a ComposedExpressions structure.","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"The function pflow follows the same API, as any regular AutoDiff.AbstractExpression.","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"julia> n_balance = length(pflow)\n14\n\njulia> round.(pflow(stack); digits=6) # evaluate the power flow balance\n14-element Vector{Float64}:\n -1.63\n -0.85\n 0.0\n 0.9\n 0.0\n 1.0\n 0.0\n 1.25\n -0.167\n 0.042\n -0.2835\n 0.171\n -0.2275\n 0.259\n","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"When we evaluate a ComposedExpressions, ExaPF first computes the basis psi(v theta) inside stack.psi, and then ExaPF uses the values in stack.psi to evaluate the final result.","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"The procedure remains the same if one wants to evaluate the LineFlows or the PowerGenerationBounds. For instance, evaluating the line flows amounts to","category":"page"},{"location":"man/formulations/","page":"Polar formulation","title":"Polar formulation","text":"julia> line_flows = ExaPF.LineFlows(polar) ∘ basis;\n\njulia> line_flows(stack)\n18-element Vector{Float64}:\n 0.0\n 0.006241000000000099\n 0.0320410000000001\n 0.0\n 0.010920249999999961\n 0.005550250000000068\n 0.0\n 0.02340899999999987\n 0.007743999999999858\n 0.0\n 0.006241000000000099\n 0.0320410000000001\n 0.0\n 0.010920249999999961\n 0.005550250000000068\n 0.0\n 0.02340899999999987\n 0.007743999999999858\n","category":"page"},{"location":"quickstart/#Quick-Start","page":"Quick start","title":"Quick Start","text":"","category":"section"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"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.","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"We start by importing CUDA and KernelAbstractions:","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"using CUDA\nusing KernelAbstractions","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"Then, we load ExaPF and its submodules with","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"using ExaPF\nimport ExaPF: AutoDiff\nconst PS = ExaPF.PowerSystem\nconst LS = ExaPF.LinearSolvers","category":"page"},{"location":"quickstart/#Short-version","page":"Quick start","title":"Short version","text":"","category":"section"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"ExaPF loads instances from the pglib-opf benchmark. ExaPF contains an artifact defined in Artifacts.toml that is built from the ExaData repository containing Exascale Computing Project relevant test cases. You may set a data file using","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"datafile = joinpath(artifact\"ExaData\", \"ExaData\", \"case1354.m\")","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"using LazyArtifacts\nusing ExaPF\nusing CUDA\nusing KernelAbstractions\nusing ExaPF\nimport ExaPF: AutoDiff\nconst PS = ExaPF.PowerSystem\nconst LS = ExaPF.LinearSolvers\nartifact_toml = joinpath(@__DIR__, \"..\", \"..\", \"Artifacts.toml\")\nexadata_hash = artifact_hash(\"ExaData\", artifact_toml)\ndatafile = joinpath(artifact_path(exadata_hash), \"ExaData\", \"case1354.m\")","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"The powerflow equations can be solved in three lines of code, as","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"polar = ExaPF.PolarForm(datafile, CPU()) # Load data\nstack = ExaPF.NetworkStack(polar) # Load variables\nconvergence = run_pf(polar, stack; verbose=1)\n","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"Implicitly, ExaPF has just proceed to the following operations:","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"instantiate automatically a starting point x_0 from the variables stored in stack\ninstantiate the Jacobian of the powerflow equations using AutoDiff.\nsolve the powerflow equations iteratively, using a Newton-Raphson algorithm.","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"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.","category":"page"},{"location":"quickstart/#Detailed-version","page":"Quick start","title":"Detailed version","text":"","category":"section"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"In what follows, we detail step by step the detailed procedure to solve the powerflow equations.","category":"page"},{"location":"quickstart/#How-to-load-a-MATPOWER-instance-as-a-PowerNetwork-object?","page":"Quick start","title":"How to load a MATPOWER instance as a PowerNetwork object?","text":"","category":"section"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"We start by importing a MATPOWER instance to a ExaPF.PowerSystem.PowerNetwork object:","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"pf = PS.PowerNetwork(datafile)","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"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","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"nbus = PS.get(pf, PS.NumberOfBuses())\npv_indexes = pf.pv;","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"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.","category":"page"},{"location":"quickstart/#How-to-solve-the-powerflow-equations?","page":"Quick start","title":"How to solve the powerflow equations?","text":"","category":"section"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"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 bmv = ve^j theta, each bus i=1 cdots N_B must satisfy the power balance equations:","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"beginaligned\n p_i = v_i sum_j^n v_j (g_ijcos(theta_i - theta_j) + b_ijsin(theta_i - theta_j)) \n q_i = v_i sum_j^n v_j (g_ijsin(theta_i - theta_j) - b_ijcos(theta_i - theta_j)) \nendaligned","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"The powerflow equations rewrite in the abstract mathematical formalism:","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"g(x u) = 0","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"For a given control u, solving the powerflow equations resumes to find a state x(u) such that g(x(u) u) = 0.","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"To this goal, ExaPF.jl implements a Newton-Raphson algorithm that allows to solve the powerflow equations in a few lines of code. We first instantiate a PolarForm object to adopt a polar formulation as a model:","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"polar = ExaPF.PolarForm(pf, CPU())","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"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().","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"The Newton-Raphson solves the equation g(x u) = 0 in an iterative fashion. The algorithm solves at each step the linear equation:","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":" x_k+1 = - (nabla_x g_k)^-1 g(x_k u)","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"Hence, the algorithm requires the following elements:","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"an initial variable x_0\na function to solve efficiently the linear system (nabla_x g_k) x_k+1 = g(x_k u)\na function to evaluate the Jacobian nabla_x g_k","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"The variable x is instantiated as:","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"stack = ExaPF.NetworkStack(polar)","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"The function g is implemented using ExaPF's custom modeler:","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"basis = ExaPF.PolarBasis(polar)\npowerflow = ExaPF.PowerFlowBalance(polar) ∘ basis","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"The Jacobian nabla_x g is evaluated automatically using forward-mode AutoDiff:","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"mapx = ExaPF.mapping(polar, State());\njx = ExaPF.Jacobian(polar, powerflow, mapx)","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"The (direct) linear solver can be instantiated directly as","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"linear_solver = LS.DirectSolver(jx.J);\n","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"Let's explain further these three objects.","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"stack is a AbstractStack storing all the variables attached to the formulation polar::PolarForm.\njx is a Jacobian structure which allows the solver to compute efficiently the Jacobian of the powerflow equations nabla_x g using AutoDiff.\nlinear_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.","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"In the AutoDiff Jacobian jx, the evaluation of the Jacobian J is stored in jx.J:","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"jac = jx.J","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"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.","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"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:","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"pf_algo = NewtonRaphson(; verbose=1, tol=1e-10)","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"Then, we can solve the powerflow equations simply with","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"convergence = ExaPF.nlsolve!(pf_algo, jx, stack; linear_solver=linear_solver)","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"Here, the algorithm solves the powerflow equations in 5 iterations. The algorithm modifies the values of stack inplace, to avoid any unnecessary memory allocations.","category":"page"},{"location":"quickstart/#How-to-deport-the-computation-on-the-GPU?","page":"Quick start","title":"How to deport the computation on the GPU?","text":"","category":"section"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"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:","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"polar_gpu = ExaPF.PolarForm(pf, CUDABackend())\n","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"polar_gpu will load all the structures it needs on the GPU, to avoid unnecessary movements between the host and the device. We can load the other structures directly on the GPU with:","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"stack_gpu = ExaPF.NetworkStack(polar_gpu)\n\nbasis_gpu = ExaPF.PolarBasis(polar_gpu)\npflow_gpu = ExaPF.PowerFlowBalance(polar_gpu) ∘ basis_gpu\njx_gpu = ExaPF.Jacobian(polar_gpu, pflow_gpu, mapx)\n\nlinear_solver = LS.DirectSolver(jx_gpu.J)","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"Then, solving the powerflow equations on the GPU directly translates as","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"convergence = ExaPF.nlsolve!(pf_algo, jx_gpu, stack_gpu; linear_solver=linear_solver)","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"Note that we get exactly the same iterations as when we solve the power flow equations on the CPU.","category":"page"},{"location":"quickstart/#How-to-solve-the-linear-system-with-BICGSTAB?","page":"Quick start","title":"How to solve the linear system with BICGSTAB?","text":"","category":"section"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"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).","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"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","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"npartitions = 8;\njac_gpu = jx_gpu.J;\nprecond = BlockJacobiPreconditioner(jac_gpu, npartitions, CUDABackend());","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"You can attach the preconditioner to an BICGSTAB algorithm simply as","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"linear_solver = ExaPF.Bicgstab(jac_gpu; P=precond);\n","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"(this will use the BICGSTAB algorithm implemented in Krylov.jl).","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"We need to update accordingly the tolerance of the Newton-Raphson algorithm (the iterative solver is less accurate than the direct solver):","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"pf_algo = NewtonRaphson(; verbose=1, tol=1e-7)","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"We reset the variables to their initial values:","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"ExaPF.init!(polar_gpu, stack_gpu)","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"Then, solving the power flow with the iterative solvers directly translates to one call to nlsolve!:","category":"page"},{"location":"quickstart/","page":"Quick start","title":"Quick start","text":"convergence = ExaPF.nlsolve!(pf_algo, jx_gpu, stack_gpu; linear_solver=linear_solver)","category":"page"},{"location":"tutorials/batch_evaluation/","page":"Power flow: batch evaluation","title":"Power flow: batch evaluation","text":"CurrentModule = ExaPF\nDocTestSetup = quote\n using ExaPF\n const LS = ExaPF.LinearSolvers\n const AD = ExaPF.AD\nend\nDocTestFilters = [r\"ExaPF\"]","category":"page"},{"location":"tutorials/batch_evaluation/","page":"Power flow: batch evaluation","title":"Power flow: batch evaluation","text":"using ExaPF\nusing KLU\nusing LinearAlgebra\nconst LS = ExaPF.LinearSolvers\nconst PS = ExaPF.PowerSystem\n\npolar = ExaPF.load_polar(\"case9.m\")","category":"page"},{"location":"tutorials/batch_evaluation/#Batch-power-flow","page":"Power flow: batch evaluation","title":"Batch power flow","text":"","category":"section"},{"location":"tutorials/batch_evaluation/","page":"Power flow: batch evaluation","title":"Power flow: batch evaluation","text":"ExaPF provides a way to evaluate the expressions by blocks, opening the way to introduce more parallelism in the code.","category":"page"},{"location":"tutorials/batch_evaluation/#BlockPolarForm","page":"Power flow: batch evaluation","title":"BlockPolarForm","text":"","category":"section"},{"location":"tutorials/batch_evaluation/","page":"Power flow: batch evaluation","title":"Power flow: batch evaluation","text":"We recall that a given NetworkStack stack stores the different variables and parameters (power generations, voltages, loads) required to evaluate the power flow model.","category":"page"},{"location":"tutorials/batch_evaluation/","page":"Power flow: batch evaluation","title":"Power flow: batch evaluation","text":"stack = ExaPF.NetworkStack(polar);\n","category":"page"},{"location":"tutorials/batch_evaluation/","page":"Power flow: batch evaluation","title":"Power flow: batch evaluation","text":"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","category":"page"},{"location":"tutorials/batch_evaluation/","page":"Power flow: batch evaluation","title":"Power flow: batch evaluation","text":"nbus = ExaPF.get(polar, PS.NumberOfBuses());\npd = stack.params[1:nbus]\nqd = stack.params[nbus+1:2*nbus]\n","category":"page"},{"location":"tutorials/batch_evaluation/","page":"Power flow: batch evaluation","title":"Power flow: batch evaluation","text":"By default, a NetworkStack stores one set of loads p_0.","category":"page"},{"location":"tutorials/batch_evaluation/","page":"Power flow: batch evaluation","title":"Power flow: batch evaluation","text":"Suppose now we want to evaluate the model associated with the polar formulation for N different set of parameters (=scenarios) p_1 cdots p_N. ExaPF allows to streamline the polar formulation with a BlockPolarForm structure:","category":"page"},{"location":"tutorials/batch_evaluation/","page":"Power flow: batch evaluation","title":"Power flow: batch evaluation","text":"nscen = 10;\nblk_polar = ExaPF.BlockPolarForm(polar, nscen)\n","category":"page"},{"location":"tutorials/batch_evaluation/","page":"Power flow: batch evaluation","title":"Power flow: batch evaluation","text":"Then, ExaPF can also instantiate a NetworkStack object, with the memory required to store the variables of the different scenarios:","category":"page"},{"location":"tutorials/batch_evaluation/","page":"Power flow: batch evaluation","title":"Power flow: batch evaluation","text":"blk_stack = ExaPF.NetworkStack(blk_polar)\n","category":"page"},{"location":"tutorials/batch_evaluation/","page":"Power flow: batch evaluation","title":"Power flow: batch evaluation","text":"We can pass the scenarios manually using the function set_params!:","category":"page"},{"location":"tutorials/batch_evaluation/","page":"Power flow: batch evaluation","title":"Power flow: batch evaluation","text":"\nploads = rand(nbus, nscen);\nqloads = rand(nbus, nscen);\nExaPF.set_params!(blk_stack, ploads, qloads)\n","category":"page"},{"location":"tutorials/batch_evaluation/","page":"Power flow: batch evaluation","title":"Power flow: batch evaluation","text":"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):","category":"page"},{"location":"tutorials/batch_evaluation/","page":"Power flow: batch evaluation","title":"Power flow: batch evaluation","text":"reshape(blk_stack.vmag, nbus, nscen)","category":"page"},{"location":"tutorials/batch_evaluation/","page":"Power flow: batch evaluation","title":"Power flow: batch evaluation","text":"Only the parameters are varying according to the scenarios we passed as input in the constructor:","category":"page"},{"location":"tutorials/batch_evaluation/","page":"Power flow: batch evaluation","title":"Power flow: batch evaluation","text":"reshape(blk_stack.pload, nbus, nscen)","category":"page"},{"location":"tutorials/batch_evaluation/#Evaluate-expressions-in-block","page":"Power flow: batch evaluation","title":"Evaluate expressions in block","text":"","category":"section"},{"location":"tutorials/batch_evaluation/","page":"Power flow: batch evaluation","title":"Power flow: batch evaluation","text":"ExaPF takes advantage of the block structure when using a BlockPolarForm.","category":"page"},{"location":"tutorials/batch_evaluation/","page":"Power flow: batch evaluation","title":"Power flow: batch evaluation","text":"As an example, suppose we want to evaluate the power flow balances in block form with a PowerFlowBalance expression:","category":"page"},{"location":"tutorials/batch_evaluation/","page":"Power flow: batch evaluation","title":"Power flow: batch evaluation","text":"powerflow = ExaPF.PowerFlowBalance(blk_polar) ∘ ExaPF.PolarBasis(blk_polar);\n","category":"page"},{"location":"tutorials/batch_evaluation/","page":"Power flow: batch evaluation","title":"Power flow: batch evaluation","text":"A block evaluation takes as input the NetworkStack blk_stack structure:","category":"page"},{"location":"tutorials/batch_evaluation/","page":"Power flow: batch evaluation","title":"Power flow: batch evaluation","text":"m = div(length(powerflow), nscen);\nblk_output = zeros(m * nscen);\npowerflow(blk_output, blk_stack); # inplace evaluation\nreshape(blk_output, m, nscen)\n","category":"page"},{"location":"tutorials/batch_evaluation/","page":"Power flow: batch evaluation","title":"Power flow: batch evaluation","text":"We get N different results for the power flow balance equations, depending on which scenario we are on.","category":"page"},{"location":"tutorials/batch_evaluation/#Solve-power-flow-in-block-on-the-CPU","page":"Power flow: batch evaluation","title":"Solve power flow in block on the CPU","text":"","category":"section"},{"location":"tutorials/batch_evaluation/","page":"Power flow: batch evaluation","title":"Power flow: batch evaluation","text":"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:","category":"page"},{"location":"tutorials/batch_evaluation/","page":"Power flow: batch evaluation","title":"Power flow: batch evaluation","text":"blk_jx = ExaPF.ArrowheadJacobian(blk_polar, powerflow, State());\nblk_jx.J","category":"page"},{"location":"tutorials/batch_evaluation/","page":"Power flow: batch evaluation","title":"Power flow: batch evaluation","text":"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:","category":"page"},{"location":"tutorials/batch_evaluation/","page":"Power flow: batch evaluation","title":"Power flow: batch evaluation","text":"ExaPF.set_params!(blk_jx, blk_stack);\nExaPF.jacobian!(blk_jx, blk_stack);\n","category":"page"},{"location":"tutorials/batch_evaluation/","page":"Power flow: batch evaluation","title":"Power flow: batch evaluation","text":"As soon as the blk_jx initialized, we can solve the power flow equations in block as","category":"page"},{"location":"tutorials/batch_evaluation/","page":"Power flow: batch evaluation","title":"Power flow: batch evaluation","text":"conv = ExaPF.nlsolve!(\n NewtonRaphson(verbose=2),\n blk_jx,\n blk_stack;\n)","category":"page"},{"location":"tutorials/batch_evaluation/","page":"Power flow: batch evaluation","title":"Power flow: batch evaluation","text":"At the solution, we get different values for the voltage magnitudes at the PQ nodes:","category":"page"},{"location":"tutorials/batch_evaluation/","page":"Power flow: batch evaluation","title":"Power flow: batch evaluation","text":"reshape(blk_stack.vmag, nbus, nscen)","category":"page"},{"location":"tutorials/batch_evaluation/#Solve-power-flow-in-batch-on-the-GPU","page":"Power flow: batch evaluation","title":"Solve power flow in batch on the GPU","text":"","category":"section"},{"location":"tutorials/batch_evaluation/","page":"Power flow: batch evaluation","title":"Power flow: batch evaluation","text":"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:","category":"page"},{"location":"tutorials/batch_evaluation/","page":"Power flow: batch evaluation","title":"Power flow: batch evaluation","text":"using CUDA, CUSOLVERRF\npolar_gpu = ExaPF.load_polar(\"case9.m\", CUDABackend());\nblk_polar_gpu = ExaPF.BlockPolarForm(polar_gpu, nscen); # load model on GPU\nblk_stack_gpu = ExaPF.NetworkStack(blk_polar_gpu);\nExaPF.set_params!(blk_stack_gpu, ploads, qloads);\npowerflow_gpu = ExaPF.PowerFlowBalance(blk_polar_gpu) ∘ ExaPF.PolarBasis(blk_polar_gpu);\nblk_jx_gpu = ExaPF.ArrowheadJacobian(blk_polar_gpu, powerflow_gpu, State());\nExaPF.set_params!(blk_jx_gpu, blk_stack_gpu);\nExaPF.jacobian!(blk_jx_gpu, blk_stack_gpu);\nrf_fac = CUSOLVERRF.RFLU(blk_jx_gpu.J)\nrf_solver = LS.DirectSolver(rf_fac)\nconv = ExaPF.nlsolve!(\n NewtonRaphson(verbose=2),\n blk_jx_gpu,\n blk_stack_gpu;\n linear_solver=rf_solver,\n)\n","category":"page"},{"location":"lib/powersystem/","page":"PowerSystem","title":"PowerSystem","text":"CurrentModule = ExaPF.PowerSystem","category":"page"},{"location":"lib/powersystem/#PowerSystem","page":"PowerSystem","title":"PowerSystem","text":"","category":"section"},{"location":"lib/powersystem/#Description","page":"PowerSystem","title":"Description","text":"","category":"section"},{"location":"lib/powersystem/","page":"PowerSystem","title":"PowerSystem","text":"AbstractPowerSystem\nPowerNetwork\nload_case","category":"page"},{"location":"lib/powersystem/#ExaPF.PowerSystem.AbstractPowerSystem","page":"PowerSystem","title":"ExaPF.PowerSystem.AbstractPowerSystem","text":"AbstractPowerSystem\n\nFirst layer of the package. Store the topology of a given transmission network, including:\n\nthe power injection at each bus ;\nthe admittance matrix ;\nthe default voltage at each bus.\n\nData are imported either from a matpower file, or a PSSE file.\n\n\n\n\n\n","category":"type"},{"location":"lib/powersystem/#ExaPF.PowerSystem.PowerNetwork","page":"PowerSystem","title":"ExaPF.PowerSystem.PowerNetwork","text":"PowerNetwork <: AbstractPowerSystem\n\nThis structure contains constant parameters that define the topology and physics of the power network.\n\nThe 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.\n\nNote\n\nThe object PowerNetwork is created in the host memory. Use a AbstractFormulation object to move data to the target device.\n\n\n\n\n\n","category":"type"},{"location":"lib/powersystem/#ExaPF.PowerSystem.load_case","page":"PowerSystem","title":"ExaPF.PowerSystem.load_case","text":"load_case(case_name, lib::PowerNetworkLibrary=EXADATA)\n\nConvenient 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.\n\nExamples\n\njulia> net = PS.load_case(\"case118\") # default is MATPOWER\nPowerNetwork object with:\n Buses: 118 (Slack: 1. PV: 53. PQ: 64)\n Generators: 54.\n\njulia> net = PS.load_case(\"case1354_pegase\", PS.PGLIB)\nPowerNetwork object with:\n Buses: 1354 (Slack: 1. PV: 259. PQ: 1094)\n Generators: 260.\n\n\n\n\n\n","category":"function"},{"location":"lib/powersystem/#API-Reference","page":"PowerSystem","title":"API Reference","text":"","category":"section"},{"location":"lib/powersystem/#Network-elements","page":"PowerSystem","title":"Network elements","text":"","category":"section"},{"location":"lib/powersystem/","page":"PowerSystem","title":"PowerSystem","text":"AbstractNetworkElement","category":"page"},{"location":"lib/powersystem/#ExaPF.PowerSystem.AbstractNetworkElement","page":"PowerSystem","title":"ExaPF.PowerSystem.AbstractNetworkElement","text":"AbstractNetworkElement\n\nAbstraction for all physical elements being parts of a AbstractPowerSystem. Elements are divided in\n\ntransmission lines (Lines)\nbuses (Buses)\ngenerators (Generators)\n\n\n\n\n\n","category":"type"},{"location":"lib/powersystem/","page":"PowerSystem","title":"PowerSystem","text":"List of elements:","category":"page"},{"location":"lib/powersystem/","page":"PowerSystem","title":"PowerSystem","text":"Buses\nLines\nGenerators","category":"page"},{"location":"lib/powersystem/#ExaPF.PowerSystem.Buses","page":"PowerSystem","title":"ExaPF.PowerSystem.Buses","text":"Buses <: AbstractNetworkElement\n\nBuses of a transmission network.\n\n\n\n\n\n","category":"type"},{"location":"lib/powersystem/#ExaPF.PowerSystem.Lines","page":"PowerSystem","title":"ExaPF.PowerSystem.Lines","text":"Lines <: AbstractNetworkElement\n\nLines of a transmission network.\n\n\n\n\n\n","category":"type"},{"location":"lib/powersystem/#ExaPF.PowerSystem.Generators","page":"PowerSystem","title":"ExaPF.PowerSystem.Generators","text":"Generators <: AbstractElement\n\nGenerators in a transmission network\n\n\n\n\n\n","category":"type"},{"location":"lib/powersystem/#Network-attributes","page":"PowerSystem","title":"Network attributes","text":"","category":"section"},{"location":"lib/powersystem/","page":"PowerSystem","title":"PowerSystem","text":"AbstractNetworkAttribute","category":"page"},{"location":"lib/powersystem/#ExaPF.PowerSystem.AbstractNetworkAttribute","page":"PowerSystem","title":"ExaPF.PowerSystem.AbstractNetworkAttribute","text":"AbstractNetworkAttribute\n\nAttribute of a AbstractPowerSystem.\n\n\n\n\n\n","category":"type"},{"location":"lib/powersystem/","page":"PowerSystem","title":"PowerSystem","text":"List of attributes:","category":"page"},{"location":"lib/powersystem/","page":"PowerSystem","title":"PowerSystem","text":"NumberOfBuses\nNumberOfLines\nNumberOfGenerators\nNumberOfPVBuses\nNumberOfPQBuses\nNumberOfSlackBuses\nBaseMVA\nBusAdmittanceMatrix","category":"page"},{"location":"lib/powersystem/#ExaPF.PowerSystem.NumberOfBuses","page":"PowerSystem","title":"ExaPF.PowerSystem.NumberOfBuses","text":"NumberOfBuses <: AbstractNetworkAttribute\n\nNumber of buses in a AbstractPowerSystem.\n\n\n\n\n\n","category":"type"},{"location":"lib/powersystem/#ExaPF.PowerSystem.NumberOfLines","page":"PowerSystem","title":"ExaPF.PowerSystem.NumberOfLines","text":"NumberOfLines <: AbstractNetworkAttribute\n\nNumber of lines in a AbstractPowerSystem.\n\n\n\n\n\n","category":"type"},{"location":"lib/powersystem/#ExaPF.PowerSystem.NumberOfGenerators","page":"PowerSystem","title":"ExaPF.PowerSystem.NumberOfGenerators","text":"NumberOfGenerators <: AbstractNetworkAttribute\n\nNumber of generators in a AbstractPowerSystem.\n\n\n\n\n\n","category":"type"},{"location":"lib/powersystem/#ExaPF.PowerSystem.NumberOfPVBuses","page":"PowerSystem","title":"ExaPF.PowerSystem.NumberOfPVBuses","text":"NumberOfPVBuses <: AbstractNetworkAttribute\n\nNumber of PV buses in a AbstractPowerSystem.\n\n\n\n\n\n","category":"type"},{"location":"lib/powersystem/#ExaPF.PowerSystem.NumberOfPQBuses","page":"PowerSystem","title":"ExaPF.PowerSystem.NumberOfPQBuses","text":"NumberOfPQBuses <: AbstractNetworkAttribute\n\nNumber of PQ buses in a AbstractPowerSystem.\n\n\n\n\n\n","category":"type"},{"location":"lib/powersystem/#ExaPF.PowerSystem.NumberOfSlackBuses","page":"PowerSystem","title":"ExaPF.PowerSystem.NumberOfSlackBuses","text":"NumberOfSlackBuses <: AbstractNetworkAttribute\n\nNumber of slack buses in a AbstractPowerSystem.\n\n\n\n\n\n","category":"type"},{"location":"lib/powersystem/#ExaPF.PowerSystem.BaseMVA","page":"PowerSystem","title":"ExaPF.PowerSystem.BaseMVA","text":"BaseMVA <: AbstractNetworkAttribute\n\nBase MVA of the network.\n\n\n\n\n\n","category":"type"},{"location":"lib/powersystem/#ExaPF.PowerSystem.BusAdmittanceMatrix","page":"PowerSystem","title":"ExaPF.PowerSystem.BusAdmittanceMatrix","text":"BusAdmittanceMatrix <: AbstractNetworkAttribute\n\nBus admittance matrix associated with the topology of the network.\n\n\n\n\n\n","category":"type"},{"location":"lib/powersystem/","page":"PowerSystem","title":"PowerSystem","text":"Query the indexing of the different elements in a given network:","category":"page"},{"location":"lib/powersystem/","page":"PowerSystem","title":"PowerSystem","text":"PVIndexes\nPQIndexes\nSlackIndexes\nGeneratorIndexes\n","category":"page"},{"location":"lib/powersystem/#ExaPF.PowerSystem.PVIndexes","page":"PowerSystem","title":"ExaPF.PowerSystem.PVIndexes","text":"PVIndexes <: AbstractIndexing\n\nIndexes of the PV buses in a AbstractPowerSystem.\n\n\n\n\n\n","category":"type"},{"location":"lib/powersystem/#ExaPF.PowerSystem.PQIndexes","page":"PowerSystem","title":"ExaPF.PowerSystem.PQIndexes","text":"PQIndexes <: AbstractIndexing\n\nIndexes of the PQ buses in a AbstractPowerSystem.\n\n\n\n\n\n","category":"type"},{"location":"lib/powersystem/#ExaPF.PowerSystem.SlackIndexes","page":"PowerSystem","title":"ExaPF.PowerSystem.SlackIndexes","text":"SlackIndexes <: AbstractIndexing\n\nIndexes of the slack buses in a AbstractPowerSystem.\n\n\n\n\n\n","category":"type"},{"location":"lib/powersystem/#ExaPF.PowerSystem.GeneratorIndexes","page":"PowerSystem","title":"ExaPF.PowerSystem.GeneratorIndexes","text":"GeneratorIndexes <: AbstractIndexing\n\nIndexes of the generators in a AbstractPowerSystem.\n\n\n\n\n\n","category":"type"},{"location":"lib/powersystem/#Network-values","page":"PowerSystem","title":"Network values","text":"","category":"section"},{"location":"lib/powersystem/","page":"PowerSystem","title":"PowerSystem","text":"AbstractNetworkValues","category":"page"},{"location":"lib/powersystem/#ExaPF.PowerSystem.AbstractNetworkValues","page":"PowerSystem","title":"ExaPF.PowerSystem.AbstractNetworkValues","text":"AbstractNetworkValues\n\nNumerical values attached to the different attributes of the network.\n\n\n\n\n\n","category":"type"},{"location":"lib/powersystem/","page":"PowerSystem","title":"PowerSystem","text":"List of values:","category":"page"},{"location":"lib/powersystem/","page":"PowerSystem","title":"PowerSystem","text":"VoltageMagnitude\nVoltageAngle\nActivePower\nReactivePower\nActiveLoad\nReactiveLoad\n","category":"page"},{"location":"lib/powersystem/#ExaPF.PowerSystem.VoltageMagnitude","page":"PowerSystem","title":"ExaPF.PowerSystem.VoltageMagnitude","text":"VoltageMagnitude <: AbstractNetworkValues\n\nMagnitude |v| of the voltage v = |v| exp(i θ).\n\n\n\n\n\n","category":"type"},{"location":"lib/powersystem/#ExaPF.PowerSystem.VoltageAngle","page":"PowerSystem","title":"ExaPF.PowerSystem.VoltageAngle","text":"VoltageAngle <: AbstractNetworkValues\n\nAngle θ of the voltage v = |v| exp(i θ).\n\n\n\n\n\n","category":"type"},{"location":"lib/powersystem/#ExaPF.PowerSystem.ActivePower","page":"PowerSystem","title":"ExaPF.PowerSystem.ActivePower","text":"ActivePower <: AbstractNetworkValues\n\nActive power P of the complex power S = P + iQ.\n\n\n\n\n\n","category":"type"},{"location":"lib/powersystem/#ExaPF.PowerSystem.ReactivePower","page":"PowerSystem","title":"ExaPF.PowerSystem.ReactivePower","text":"ReactivePower <: AbstractNetworkValues\n\nReactive power Q of the complex power S = P + iQ.\n\n\n\n\n\n","category":"type"},{"location":"lib/powersystem/#ExaPF.PowerSystem.ActiveLoad","page":"PowerSystem","title":"ExaPF.PowerSystem.ActiveLoad","text":"ActiveLoad <: AbstractNetworkValues\n\nActive load Pd at buses.\n\n\n\n\n\n","category":"type"},{"location":"lib/powersystem/#ExaPF.PowerSystem.ReactiveLoad","page":"PowerSystem","title":"ExaPF.PowerSystem.ReactiveLoad","text":"ReactiveLoad <: AbstractNetworkValues\n\nReactive load Qd at buses.\n\n\n\n\n\n","category":"type"},{"location":"lib/powersystem/","page":"PowerSystem","title":"PowerSystem","text":"Function to get the range of a given value:","category":"page"},{"location":"lib/powersystem/","page":"PowerSystem","title":"PowerSystem","text":"bounds","category":"page"},{"location":"lib/powersystem/#ExaPF.PowerSystem.bounds","page":"PowerSystem","title":"ExaPF.PowerSystem.bounds","text":"bounds(pf::AbstractPowerSystem, attr::AbstractNetworkAttribute, val::AbstractNetworkValues)\n\nReturn lower and upper bounds corresponding to the admissible values of the AbstractNetworkAttribute attr.\n\nExamples\n\np_min, p_max = bounds(pf, Generator(), ActivePower())\nv_min, v_max = bounds(pf, Buses(), VoltageMagnitude())\n\n\n\n\n\n\n","category":"function"},{"location":"#ExaPF","page":"Home","title":"ExaPF","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"ExaPF.jl is a package to solve the power flow problem on upcoming exascale architectures. The code has been designed to be:","category":"page"},{"location":"","page":"Home","title":"Home","text":"Portable: Targeting exascale architectures implies a focus on graphics processing units (GPUs) as these systems lack substantial computational performance through classical CPUs.\nDifferentiable: All the expressions implemented in ExaPF are fully compatible with ForwardDiff.jl, and routines are provided to extract first- and second-order derivatives to solve efficiently power flow and optimal power flow problems.","category":"page"},{"location":"","page":"Home","title":"Home","text":"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 provide extensions that leverage the packages CUDA.jl, [AMDGPU.jl]((https://github.com/JuliaGPU/AMDGPU.jl), and KernelAbstractions.jl to make ExaPF portable across GPU architectures.","category":"page"},{"location":"#Table-of-contents","page":"Home","title":"Table of contents","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"Pages = [\n \"quickstart.md\",\n]\nDepth=1","category":"page"},{"location":"#Manual","page":"Home","title":"Manual","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"Pages = [\n \"man/formulations.md\",\n \"man/powersystem.md\",\n \"man/autodiff.md\",\n \"man/linearsolver.md\",\n \"man/benchmark.md\",\n]\nDepth = 1","category":"page"},{"location":"#Library","page":"Home","title":"Library","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"Pages = [\n \"lib/formulations.md\",\n \"lib/powersystem.md\",\n \"lib/autodiff.md\",\n \"lib/linearsolver.md\",\n]\nDepth = 1","category":"page"},{"location":"#Artifact","page":"Home","title":"Artifact","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"Pages = [\n \"artifact.md\",\n]\nDepth = 1","category":"page"},{"location":"#Funding","page":"Home","title":"Funding","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"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.","category":"page"},{"location":"lib/linearsolver/","page":"Linear Solvers","title":"Linear Solvers","text":"CurrentModule = ExaPF.LinearSolvers","category":"page"},{"location":"lib/linearsolver/#Linear-solvers","page":"Linear Solvers","title":"Linear solvers","text":"","category":"section"},{"location":"lib/linearsolver/#Description","page":"Linear Solvers","title":"Description","text":"","category":"section"},{"location":"lib/linearsolver/","page":"Linear Solvers","title":"Linear Solvers","text":"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!.","category":"page"},{"location":"lib/linearsolver/","page":"Linear Solvers","title":"Linear Solvers","text":"ldiv!","category":"page"},{"location":"lib/linearsolver/#ExaPF.LinearSolvers.ldiv!","page":"Linear Solvers","title":"ExaPF.LinearSolvers.ldiv!","text":"ldiv!(solver, x, A, y)\nldiv!(solver, x, y)\n\nsolver::AbstractLinearSolver: linear solver to solve the system\nx::AbstractVector: Solution\nA::AbstractMatrix: Input matrix\ny::AbstractVector: RHS\n\nSolve 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.\n\n\n\n\n\n","category":"function"},{"location":"lib/linearsolver/#Direct-solvers","page":"Linear Solvers","title":"Direct solvers","text":"","category":"section"},{"location":"lib/linearsolver/","page":"Linear Solvers","title":"Linear Solvers","text":"ExaPF wraps UMFPACK (shipped with SuiteSparse.jl) on the CPU, and CUSPARSE on CUDA device.","category":"page"},{"location":"lib/linearsolver/","page":"Linear Solvers","title":"Linear Solvers","text":"DirectSolver","category":"page"},{"location":"lib/linearsolver/#ExaPF.LinearSolvers.DirectSolver","page":"Linear Solvers","title":"ExaPF.LinearSolvers.DirectSolver","text":"DirectSolver <: AbstractLinearSolver\n\nSolve linear system A x = y with direct linear algebra.\n\nOn the CPU, DirectSolver uses UMFPACK to solve the linear system\nOn CUDA GPU, DirectSolver redirects the resolution to the method CUSOLVER.csrlsvqr\n\n\n\n\n\n","category":"type"},{"location":"lib/linearsolver/#Iterative-solvers","page":"Linear Solvers","title":"Iterative solvers","text":"","category":"section"},{"location":"lib/linearsolver/","page":"Linear Solvers","title":"Linear Solvers","text":"Bicgstab\nDqgmres","category":"page"},{"location":"lib/linearsolver/#ExaPF.LinearSolvers.Bicgstab","page":"Linear Solvers","title":"ExaPF.LinearSolvers.Bicgstab","text":"Bicgstab <: AbstractIterativeLinearSolver\nBicgstab(precond; verbose=0, rtol=1e-10, atol=1e-10)\n\nWrap Krylov.jl's BICGSTAB algorithm to solve iteratively the linear system A x = y.\n\n\n\n\n\n","category":"type"},{"location":"lib/linearsolver/#ExaPF.LinearSolvers.Dqgmres","page":"Linear Solvers","title":"ExaPF.LinearSolvers.Dqgmres","text":"Dqgmres <: AbstractIterativeLinearSolver\nDqgmres(precond; verbose=false, memory=4)\n\nWrap Krylov.jl's Dqgmres algorithm to solve iteratively the linear system A x = y.\n\n\n\n\n\n","category":"type"},{"location":"lib/linearsolver/","page":"Linear Solvers","title":"Linear Solvers","text":"Available linear solvers can be queried with","category":"page"},{"location":"lib/linearsolver/","page":"Linear Solvers","title":"Linear Solvers","text":"list_solvers\n","category":"page"},{"location":"lib/linearsolver/#ExaPF.LinearSolvers.list_solvers","page":"Linear Solvers","title":"ExaPF.LinearSolvers.list_solvers","text":"list_solvers(::KernelAbstractions.Device)\n\nList linear solvers available on current device.\n\n\n\n\n\n","category":"function"},{"location":"lib/linearsolver/","page":"Linear Solvers","title":"Linear Solvers","text":"A default solver is provided for each vendor backend.","category":"page"},{"location":"lib/linearsolver/","page":"Linear Solvers","title":"Linear Solvers","text":"default_linear_solver\n","category":"page"},{"location":"lib/linearsolver/#ExaPF.LinearSolvers.default_linear_solver","page":"Linear Solvers","title":"ExaPF.LinearSolvers.default_linear_solver","text":"default_linear_solver(A, ::KA.CPU)\n\nDefault linear solver on the CPU.\n\n\n\n\n\n","category":"function"},{"location":"lib/linearsolver/#Block-Krylov-solvers","page":"Linear Solvers","title":"Block-Krylov solvers","text":"","category":"section"},{"location":"lib/linearsolver/","page":"Linear Solvers","title":"Linear Solvers","text":"ExaPF.jl provides an implementation of block-GMRES to solve linear systems with multiple righ-hand sides.","category":"page"},{"location":"lib/linearsolver/","page":"Linear Solvers","title":"Linear Solvers","text":"BlockKrylovSolver\nBlockGmresSolver\nblock_gmres\nblock_gmres!","category":"page"},{"location":"lib/linearsolver/#ExaPF.LinearSolvers.BlockKrylovSolver","page":"Linear Solvers","title":"ExaPF.LinearSolvers.BlockKrylovSolver","text":"Abstract type for using block Krylov solvers in-place\n\n\n\n\n\n","category":"type"},{"location":"lib/linearsolver/#ExaPF.LinearSolvers.BlockGmresSolver","page":"Linear Solvers","title":"ExaPF.LinearSolvers.BlockGmresSolver","text":"Type for storing the vectors required by the in-place version of BLOCK-GMRES.\n\nThe outer constructors\n\nsolver = BlockGmresSolver(m, n, p, memory, SV, SM)\nsolver = BlockGmresSolver(A, B; memory=5)\n\nmay be used in order to create these vectors. memory is set to div(n,p) if the value given is larger than div(n,p).\n\n\n\n\n\n","category":"type"},{"location":"lib/linearsolver/#ExaPF.LinearSolvers.block_gmres","page":"Linear Solvers","title":"ExaPF.LinearSolvers.block_gmres","text":"(X, stats) = block_gmres(A, B; X0::AbstractMatrix{FC}; memory::Int=20, M=I, N=I,\n ldiv::Bool=false, restart::Bool=false, reorthogonalization::Bool=false,\n atol::T = √eps(T), rtol::T=√eps(T), itmax::Int=0,\n timemax::Float64=Inf, verbose::Int=0, history::Bool=false)\n\nT is an AbstractFloat such as Float32, Float64 or BigFloat. FC is T or Complex{T}.\n\n(X, stats) = block_gmres(A, B, X0::AbstractVector; kwargs...)\n\nGMRES can be warm-started from an initial guess X0 where kwargs are the same keyword arguments as above.\n\nSolve the linear system AX = B of size n with p right-hand sides using block-GMRES.\n\nInput arguments\n\nA: a linear operator that models a matrix of dimension n;\nB: a matrix of size n × p.\n\nOptional argument\n\nX0: a matrix of size n × p that represents an initial guess of the solution X.\n\nKeyword arguments\n\nmemory: if restart = true, the restarted version block-GMRES(k) is used with k = memory. If restart = false, the parameter memory should be used as a hint of the number of iterations to limit dynamic memory allocations. Additional storage will be allocated if the number of iterations exceeds memory;\nM: linear operator that models a nonsingular matrix of size n used for left preconditioning;\nN: linear operator that models a nonsingular matrix of size n used for right preconditioning;\nldiv: define whether the preconditioners use ldiv! or mul!;\nrestart: restart the method after memory iterations;\nreorthogonalization: reorthogonalize the new matrices of the block-Krylov basis against all previous matrix;\natol: absolute stopping tolerance based on the residual norm;\nrtol: relative stopping tolerance based on the residual norm;\nitmax: the maximum number of iterations. If itmax=0, the default number of iterations is set to 2 * div(n,p);\ntimemax: the time limit in seconds;\nverbose: additional details can be displayed if verbose mode is enabled (verbose > 0). Information will be displayed every verbose iterations;\nhistory: collect additional statistics on the run such as residual norms.\n\nOutput arguments\n\nx: a dense matrix of size n × p;\nstats: statistics collected on the run in a BlockGmresStats.\n\n\n\n\n\n","category":"function"},{"location":"lib/linearsolver/#ExaPF.LinearSolvers.block_gmres!","page":"Linear Solvers","title":"ExaPF.LinearSolvers.block_gmres!","text":"solver = block_gmres!(solver::BlockGmresSolver, B; kwargs...)\nsolver = block_gmres!(solver::BlockGmresSolver, B, X0; kwargs...)\n\nwhere kwargs are keyword arguments of block_gmres.\n\nSee BlockGmresSolver for more details about the solver.\n\n\n\n\n\n","category":"function"},{"location":"man/autodiff/","page":"AutoDiff","title":"AutoDiff","text":"DocTestSetup = quote\n using ExaPF\n const AD = ExaPF.AD\nend\nDocTestFilters = [r\"ExaPF\"]","category":"page"},{"location":"man/autodiff/#AutoDiff","page":"AutoDiff","title":"AutoDiff","text":"","category":"section"},{"location":"man/autodiff/#Overview","page":"AutoDiff","title":"Overview","text":"","category":"section"},{"location":"man/autodiff/","page":"AutoDiff","title":"AutoDiff","text":"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.","category":"page"},{"location":"man/autodiff/","page":"AutoDiff","title":"AutoDiff","text":" go = true\n while(go)\n dx .= jacobian(x)\\f(x)\n x .= x .- dx\n go = norm(f(x)) < tol ? true : false\n end","category":"page"},{"location":"man/autodiff/","page":"AutoDiff","title":"AutoDiff","text":"There are two modes of differentiation called forward/tangent or reverse/adjoint. The latter is known in machine learning as backpropagation. The forward mode generates Jacobian-vector product code tgt(x,d) = J(x) * d, while the adjoint mode generates code for the transposed Jacobian-vector product adj(x,y) = (J(x)'*y). We recommend the book Evaluating derivatives: principles and techniques of algorithmic differentiation by Griewank and Walther[1] for a more in-depth introduction to automatic differentiation. The computational complexity of both models favors the adjoint mode if the number of outputs of F is much smaller than the number of inputs size(x) >> size(F), like for example the loss functions in machine learning. However, in our case F is a multivariate vector function from mathbbR^n to mathbbR^n, where n is the number of buses.","category":"page"},{"location":"man/autodiff/","page":"AutoDiff","title":"AutoDiff","text":"(Image: Jacobian coloring \\label{fig:coloring})","category":"page"},{"location":"man/autodiff/","page":"AutoDiff","title":"AutoDiff","text":"To avoid a complexity of mathcalO(n) cdot cost(F) by letting the tangent mode run over all Cartesian basis vectors of mathbbR^n, we apply the technique of Jacobian coloring to compress the sparse Jacobian J. Running the tangent mode, it allows to compute columns of the Jacobian concurrently, by combining independent columns in one Jacobian-vector evaluation (see in figure above). For sparsity detection we rely on the greedy algorithm implemented by SparseDiffTools.jl.","category":"page"},{"location":"man/autodiff/","page":"AutoDiff","title":"AutoDiff","text":"Given the sparsity pattern, the forward model is applied through the package ForwardDiff.jl. Given the number of Jacobian colors c we can build our dual type t1s with c directions:","category":"page"},{"location":"man/autodiff/","page":"AutoDiff","title":"AutoDiff","text":"t1s{N} = ForwardDiff.Dual{Nothing,Float64, N} where N}","category":"page"},{"location":"man/autodiff/","page":"AutoDiff","title":"AutoDiff","text":"Note that a second-order type t2s can be created naturally by applying the same logic to t1s:","category":"page"},{"location":"man/autodiff/","page":"AutoDiff","title":"AutoDiff","text":"t2s{M,N} = ForwardDiff.Dual{Nothing,t1s{N}, M} where M, N}","category":"page"},{"location":"man/autodiff/","page":"AutoDiff","title":"AutoDiff","text":"Finally, this dual type can be ported to both vector types Vector and CuVector:","category":"page"},{"location":"man/autodiff/","page":"AutoDiff","title":"AutoDiff","text":"VT = Vector{Float64}\nVT = Vector{t1s{N}}}\nVT = CuVector{t1s{N}}}","category":"page"},{"location":"man/autodiff/","page":"AutoDiff","title":"AutoDiff","text":"Setting VT to either of the three types allows us to instantiate code that has been written using the broadcast operator .","category":"page"},{"location":"man/autodiff/","page":"AutoDiff","title":"AutoDiff","text":"x .= a .* b","category":"page"},{"location":"man/autodiff/","page":"AutoDiff","title":"AutoDiff","text":"or accessed in kernels written for KernelAbstractions.jl like for example the power flow equations (here in polar form):","category":"page"},{"location":"man/autodiff/","page":"AutoDiff","title":"AutoDiff","text":"@kernel function residual_kernel!(F, v_m, v_a,\n ybus_re_nzval, ybus_re_colptr, ybus_re_rowval,\n ybus_im_nzval, ybus_im_colptr, ybus_im_rowval,\n pinj, qinj, pv, pq, nbus)\n\n npv = size(pv, 1)\n npq = size(pq, 1)\n\n i = @index(Global, Linear)\n # REAL PV: 1:npv\n # REAL PQ: (npv+1:npv+npq)\n # IMAG PQ: (npv+npq+1:npv+2npq)\n fr = (i <= npv) ? pv[i] : pq[i - npv]\n F[i] -= pinj[fr]\n if i > npv\n F[i + npq] -= qinj[fr]\n end\n @inbounds for c in ybus_re_colptr[fr]:ybus_re_colptr[fr+1]-1\n to = ybus_re_rowval[c]\n aij = v_a[fr] - v_a[to]\n coef_cos = v_m[fr]*v_m[to]*ybus_re_nzval[c]\n coef_sin = v_m[fr]*v_m[to]*ybus_im_nzval[c]\n cos_val = cos(aij)\n sin_val = sin(aij)\n F[i] += coef_cos * cos_val + coef_sin * sin_val\n if i > npv\n F[npq + i] += coef_cos * sin_val - coef_sin * cos_val\n end\n end\nend","category":"page"},{"location":"man/autodiff/","page":"AutoDiff","title":"AutoDiff","text":"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","category":"page"},{"location":"man/autodiff/","page":"AutoDiff","title":"AutoDiff","text":"(Image: SIMD AD for point-wise vector product \\label{fig:simd})","category":"page"},{"location":"man/autodiff/","page":"AutoDiff","title":"AutoDiff","text":"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.","category":"page"},{"location":"man/autodiff/","page":"AutoDiff","title":"AutoDiff","text":"[1]: Griewank, Andreas, and Andrea Walther. Evaluating derivatives: principles and techniques of algorithmic differentiation. Society for Industrial and Applied Mathematics, 2008.","category":"page"}] } diff --git a/previews/PR284/tutorials/batch_evaluation/index.html b/previews/PR284/tutorials/batch_evaluation/index.html index 073c0460..4000ea09 100644 --- a/previews/PR284/tutorials/batch_evaluation/index.html +++ b/previews/PR284/tutorials/batch_evaluation/index.html @@ -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.0336311224676642
- 0.5859502535737532
- 0.6145145416955379
- 0.002341794151535792
- 0.49396910291522844
- 0.9023901501554883
- 0.89031541632394
- 0.3575025449412078
- 0.9482064379261469
- 0.7667973128395988
+ 0.41502230349088065
+ 0.396944132428947
+ 0.5360256637644711
+ 0.007299211454480203
+ 0.04940031022456892
+ 0.23101254407798322
+ 0.39799009911986294
+ 0.5193993117885212
+ 0.1682183267400389
+ 0.7030264554549582
  ⋮
- 0.01531753770796096
- 0.3805120481631574
- 0.5661522366138103
- 0.5482677632413161
- 0.8077794653957094
- 0.9290792870731491
- 0.54836016336943
- 0.2518202628136732
- 0.19223818421785555

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.19062540829689023
+ 0.41746029394604
+ 0.16243613587786832
+ 0.8039090943873148
+ 0.8865657407617652
+ 0.45567865831476084
+ 0.1417572333683096
+ 0.6641876130005634
+ 0.006093086782478441

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,32 +50,32 @@
  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.704841  0.83734    0.639888   0.274402  …  0.350603  0.839943   0.505134
- 0.857662  0.0776943  0.0881648  0.836244     0.428564  0.653082   0.236166
- 0.394316  0.737298   0.824635   0.152304     0.418763  0.87285    0.836036
- 0.676135  0.854944   0.05893    0.20516      0.487551  0.088029   0.830987
- 0.854696  0.881841   0.363847   0.253802     0.714238  0.242061   0.565155
- 0.346735  0.956796   0.99845    0.664538  …  0.382553  0.583662   0.313864
- 0.664485  0.412578   0.352091   0.2571       0.145627  0.0538427  0.974443
- 0.55174   0.141624   0.710021   0.205874     0.268152  0.985651   0.00627822
- 0.579102  0.12763    0.129832   0.442655     0.134053  0.574124   0.141198

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.378325   0.348711  0.43127    0.502965  …  0.932903  0.0400788  0.330865
+ 0.975234   0.206627  0.749635   0.652573     0.489201  0.326433   0.510685
+ 0.586873   0.467072  0.762475   0.758771     0.696466  0.876382   0.270062
+ 0.330764   0.594492  0.287023   0.712396     0.799146  0.694638   0.439
+ 0.283471   0.605576  0.173379   0.56718      0.688993  0.274022   0.794709
+ 0.586467   0.638005  0.0967832  0.063373  …  0.650599  0.73444    0.905389
+ 0.834299   0.132252  0.875436   0.101527     0.240577  0.614626   0.828402
+ 0.438882   0.68637   0.856205   0.756741     0.148034  0.15367    0.278012
+ 0.0973482  0.207082  0.497734   0.935677     0.561047  0.423885   0.642365

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}:
- -0.772338  -1.55231    -1.54184    …  -1.20144    -0.976918   -1.39383
- -0.455684  -0.112702   -0.0253649     -0.431237    0.0228501  -0.0139637
-  0.676135   0.854944    0.05893        0.487551    0.088029    0.830987
-  0.854696   0.881841    0.363847       0.714238    0.242061    0.565155
-  0.346735   0.956796    0.99845        0.382553    0.583662    0.313864
-  0.664485   0.412578    0.352091   …   0.145627    0.0538427   0.974443
-  0.55174    0.141624    0.710021       0.268152    0.985651    0.00627822
-  0.579102   0.12763     0.129832       0.134053    0.574124    0.141198
- -0.164658   0.182291    0.428834       0.137061    0.171009    0.381268
-  0.235969   0.404392    0.547878      -0.219057   -0.118965    0.549779
-  0.61889   -0.209046   -0.179851   …  -0.0191358   0.397402    0.645579
-  0.711315   0.524018    0.535205       0.416454    0.795829    0.36936
-  0.130003   0.0926874   0.707446       0.224213    0.442402    0.0243203
-  0.707206   0.0168978   0.0753966      0.179711    0.614564   -0.0487618

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());
+ -0.654766   -1.42337   -0.880365   …  -1.1408     -1.30357    -1.11931
+ -0.263127   -0.382928  -0.0875245     -0.153534    0.0263817  -0.579938
+  0.330764    0.594492   0.287023       0.799146    0.694638    0.439
+  0.283471    0.605576   0.173379       0.688993    0.274022    0.794709
+  0.586467    0.638005   0.0967832      0.650599    0.73444     0.905389
+  0.834299    0.132252   0.875436   …   0.240577    0.614626    0.828402
+  0.438882    0.68637    0.856205       0.148034    0.15367     0.278012
+  0.0973482   0.207082   0.497734       0.561047    0.423885    0.642365
+ -0.159701    0.289273   0.725499       0.215777   -0.111753    0.636909
+ -0.2086     -0.151062   0.188197       0.0982724   0.0942624   0.628566
+ -0.0524875   0.360502  -0.0446838  …   0.31679     0.507153    0.172179
+  0.21899     0.327824   0.49121        0.482972   -0.127012   -0.0372428
+  0.291899    0.485098   0.0730057      0.039126    0.514546    0.436688
+ -0.0727817  -0.141673   0.315793      -0.0189263   0.263311   -0.234907

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:
 ⎡⡱⣮⡲⣞⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎤
 ⎢⣸⢮⣻⣾⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
@@ -124,20 +124,20 @@
 )
Power flow has converged: true
   * #iterations: 4
   * Time Jacobian (s) ........: 0.0002
-  * Time linear solver (s) ...: 0.0004
-     * update (s) ............: 0.0004
+  * Time linear solver (s) ...: 0.0000
+     * update (s) ............: 0.0000
      * ldiv (s) ..............: 0.0000
-  * Time total (s) ...........: 0.0008
+  * Time total (s) ...........: 0.0004

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.947387  0.963314  0.944251  0.975215     0.981789  0.953965  0.948442
- 0.909342  0.923153  0.900902  0.951059     0.9858    0.950273  0.895798
- 0.929637  0.968724  0.953182  0.984822  …  0.985741  0.939527  0.935416
- 0.893348  0.937488  0.909573  0.97681      0.957847  0.890992  0.926485
- 0.934679  0.967524  0.930956  0.991122     0.972787  0.924433  0.964219
- 0.890354  0.960474  0.932639  0.96614      0.962938  0.897604  0.953522

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
+ 1.01005   0.981236  0.937541  0.99301      0.972289  0.980127  0.939169
+ 1.01097   0.976823  0.928528  0.972866     0.951045  0.961862  0.885807
+ 0.990187  0.967132  0.963728  0.967575  …  0.960841  0.959969  0.957089
+ 0.968027  0.94775   0.928703  0.966958     0.942912  0.961512  0.955222
+ 0.977983  0.96266   0.952933  0.98836      0.973002  0.963367  0.962246
+ 1.00288   0.979431  0.916045  0.985176     0.966997  0.953981  0.953916

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);
@@ -156,8 +156,8 @@
 )
Power flow has converged: true
   * #iterations: 4
   * Time Jacobian (s) ........: 0.0014
-  * Time linear solver (s) ...: 0.0490
-     * update (s) ............: 0.0245
-     * ldiv (s) ..............: 0.0245
-  * Time total (s) ...........: 0.0524
-
+ * Time linear solver (s) ...: 0.0486 + * update (s) ............: 0.0242 + * ldiv (s) ..............: 0.0244 + * Time total (s) ...........: 0.0518 + diff --git a/previews/PR284/tutorials/direct_solver/index.html b/previews/PR284/tutorials/direct_solver/index.html index 21a6ab69..7d850f11 100644 --- a/previews/PR284/tutorials/direct_solver/index.html +++ b/previews/PR284/tutorials/direct_solver/index.html @@ -15,20 +15,20 @@ 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.1782
-  * Time linear solver (s) ...: 0.4444
-     * update (s) ............: 0.4400
-     * ldiv (s) ..............: 0.0044
-  * Time total (s) ...........: 1.6289
+  * Time Jacobian (s) ........: 1.1292
+  * Time linear solver (s) ...: 0.0337
+     * update (s) ............: 0.0294
+     * ldiv (s) ..............: 0.0043
+  * Time total (s) ...........: 1.1688

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} @0x000000000fe46670, Ptr{Nothing} @0x0000000010ed6690, 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
+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} @0x0000000011af2da0, Ptr{Nothing} @0x0000000010481480, 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.40062041224643, 995.2128749049238, -479.4136240982223, -138.29482949968133, -145.95963146384614, -230.91975549270308, 22.972728768572996  …  -480.47843944268664, -39.71137052069115, 520.1642960016878, -3360.900457400674, -208.11090120000475, 3568.982258670069, -39.40403545216686, 39.40522422897141, -242.8021538639664, 242.80196093598863]))

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.2755
-  * Time linear solver (s) ...: 0.0339
-     * update (s) ............: 0.0295
-     * ldiv (s) ..............: 0.0044
-  * Time total (s) ...........: 0.3152
+  * Time Jacobian (s) ........: 0.2735
+  * Time linear solver (s) ...: 0.0336
+     * update (s) ............: 0.0294
+     * ldiv (s) ..............: 0.0042
+  * Time total (s) ...........: 0.3131

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)
@@ -38,9 +38,9 @@
 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.6574
-  * Time linear solver (s) ...: 0.1038
-     * update (s) ............: 0.1022
-     * ldiv (s) ..............: 0.0016
-  * Time total (s) ...........: 1.7641
-
+ * Time Jacobian (s) ........: 1.6556 + * Time linear solver (s) ...: 0.1176 + * update (s) ............: 0.1161 + * ldiv (s) ..............: 0.0015 + * Time total (s) ...........: 1.7761 +