Formulations
ExaPF's formalism is based on a vectorized formulation of the power flow problem, as introduced in Lee, Turitsyn, Molzahn, Roald (2020). Throughout this page, we will refer to this formulation as LTMR2020. It is equivalent to the classical polar formulation of the OPF.
In what follows, we denote by $v \in \mathbb{R}^{n_b}$ the voltage magnitudes, $\theta \in \mathbb{R}^{n_b}$ the voltage angles and $p_g, q_g \in \mathbb{R}^{n_g}$ the active and reactive power generations. The active and reactive loads are denoted respectively by $p_d, q_d \in \mathbb{R}^{n_b}$.
Power flow model
The idea is to factorize all nonlinearities inside a basis function depending both on the voltage magnitudes $v$ and voltage angles $\theta$, such that $\psi: \mathbb{R}^{n_b} \times \mathbb{R}^{n_b} \to \mathbb{R}^{2n_\ell + n_b}$. If we introduce the intermediate expressions
\[ \psi_\ell^C(v, \theta) = v^f v^t \cos(\theta_f - \theta_t) \quad \forall \ell = 1, \cdots, n_\ell \\ \psi_\ell^S(v, \theta) = v^f v^t \sin(\theta_f - \theta_t) \quad \forall \ell = 1, \cdots, n_\ell \\ \psi_k(v, \theta) = v_k^2 \quad \forall k = 1, \cdots, n_b\]
the basis $\psi$ is defined as
\[ \psi(v, \theta) = [\psi_\ell^C(v, \theta)^\top ~ \psi_\ell^S(v, \theta)^\top ~ \psi_k(v, \theta)^\top ] \, .\]
The basis $\psi$ encodes all the nonlinearities in the problem. For instance, the power flow equations rewrite directly as
\[ \begin{bmatrix} C_g p_g - p_d \\ C_g q_g - q_d \end{bmatrix} + \underbrace{ \begin{bmatrix} - \hat{G}^c & - \hat{B}^s & -G^d \\ \hat{B}^c & - \hat{G}^s & B^d \end{bmatrix} }_{M} \psi(v, \theta) = 0\]
with $C_g \in \mathbb{R}^{n_b \times n_g}$ the bus-generators incidence matrix, and the matrices $B, G$ defined from the admittance matrix $Y_b$ of the network.
Similarly, the line flows rewrite
\[ \begin{bmatrix} s_p^f \\ s_q^f \end{bmatrix} = \overbrace{ \begin{bmatrix} G_{ft} & B_{ft} & G_{ff} C_f^\top \\ -B_{ft} & G_{ft} & -B_{ff} C_f^\top \end{bmatrix} }^{L_{line}^f} \psi(v, \theta) \\ \begin{bmatrix} s_p^t \\ s_q^t \end{bmatrix} = \underbrace{ \begin{bmatrix} G_{tf} & B_{tf} & G_{tt} C_t^\top \\ -B_{tf} & G_{tf} & -B_{tt} C_t^\top \end{bmatrix} }_{L_{line}^t} \psi(v, \theta)\]
with $C_f \in \mathbb{R}^{n_b \times n_\ell}$ the bus-from incidence matrix and $C_t \in \mathbb{R}^{n_b \times n_\ell}$ the bus-to incidence matrix. Then, the line flows constraints write directly with the quadratic expressions:
\[ (s_p^f)^2 + (s_q^f)^2 \leq (s^{max})^2 \quad \, , (s_p^t)^2 + (s_q^t)^2 \leq (s^{max})^2 \quad \, .\]
Why is this model advantageous?
Implementing the model LTMR2020 is not difficult once the basis function $\psi$ has been defined. Indeed, if we select a subset of the power flow equations (as usual, associated to the active injections at PV nodes, and active and reactive injections at PQ nodes), we get
\[ C_{eq} p_g + M_{eq} \psi + \tau = 0\]
with $C_{eq}$ defined from the bus-generator incidence matrix $C_g$, $M_{eq}$ a subset of the matrix $M$, $\tau$ a constant depending on the loads in the problem. Note that $C_{eq}$ and $M_{eq}$ are sparse matrices, so the expression can be implemented efficiently with sparse linear algebra operations (2 SpMV operations, 2 vector additions). The same holds true for the line flow constraints, evaluated with 2 SpMV operations:
\[ s^f = L_{line}^f \psi \, , \quad s^t = L_{line}^t \psi \, .\]
In ExaPF, all nonlinear expressions are written as linear operations depending on the nonlinear basis $\psi$. By doing so, all the unstructured sparsity of the power flow problem is directly handled inside the sparse linear algebra library (cusparse
on CUDA GPU, SuiteSparse
on the CPU).
In what follows, we detail the implementation of the LTMR2020 model in ExaPF.
How to instantiate the inputs?
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
.
Specifying inputs in NetworkStack
Our three inputs are $(v, \theta, p_g) \in \mathbb{R}^{2n_b + n_g}$ (voltage magnitude, voltage angle, power generations). The basis $\psi$ is considered as an intermediate expression.
We store all inputs in a NetworkStack
structure:
struct NetworkStack{VT} <: AbstractStack
input::VT
vmag::VT # voltage magnitudes (view)
vang::VT # voltage angles (view)
pgen::VT # active power generations (view)
ψ::VT # nonlinear basis ψ(vmag, vang)
end
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]
:
# Define dimension of the problem
julia> nbus, ngen, nlines = 3, 2, 4;
julia> stack = ExaPF.NetworkStack(nbus, ngen, nlines, 0, Vector{Float64}, Vector{Float64})
8-elements NetworkStack{Vector{Float64}}
julia> stack.input
8-element Vector{Float64}:
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
julia> stack.vmag .= 1;
julia> stack.vang .= 2;
julia> stack.pgen .= 3;
julia> stack.input
8-element Vector{Float64}:
1.0
1.0
1.0
2.0
2.0
2.0
3.0
3.0
The basis vector ψ
is an intermediate expression, whose values depend on the inputs.
Defining a state and a control
In the reduced space method, we have to split the variables in a state $x$ and a control $u$. By default, we define
\[ x = (\theta_{pv}, \theta_{pq}, v_{pq}) \, , \quad x = (v_{ref}, v_{pv}, p_{g,genpv}) \,.\]
and the control, and was not flexible. In the new implementation, we define the state and the control as two mappings porting on the vector stack.input
(which itself stores all the inputs in the problem):
julia> nbus, ngen, nlines = 4, 3, 4;
julia> stack = ExaPF.NetworkStack(nbus, ngen, nlines, 0, Vector{Float64}, Vector{Float64});
julia> stack.input .= 1:length(stack.input); # index array input
julia> ref, pv, pq, genpv = [1], [2], [3, 4], [2, 3];
julia> mapx = [nbus .+ pv; nbus .+ pq; pq];
julia> mapu = [ref; pv; 2*nbus .+ genpv];
julia> x = @view stack.input[mapx]
5-element view(::Vector{Float64}, [6, 7, 8, 3, 4]) with eltype Float64:
6.0
7.0
8.0
3.0
4.0
julia> u = @view stack.input[mapu]
4-element view(::Vector{Float64}, [1, 2, 10, 11]) with eltype Float64:
1.0
2.0
10.0
11.0
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.
How to manipulate the expressions?
ExaPF implements the different functions required to implement the optimal power flow problem with the polar formulation:
PowerFlowBalance
: power flow balance equationsPowerGenerationBounds
: bounds on the power generationLineFlows
: line flow constraintsCostFunction
: quadratic cost function
Each function follows the LTMR2020 model and depends on the basis function $\psi(v, \theta)$, here implemented in the PolarBasis
function.
We demonstrate how to use the different functions on the case9
instance. The procedure remains the same for all power network.
julia> polar = ExaPF.load_polar("case9.m");
julia> stack = ExaPF.NetworkStack(polar);
All the code presented below is agnostic with regards to the specific device (CPU
, CUDABackend
...) we are using. By default, ExaPF computes the expressions on the CPU. Deporting the computation on a CUDABackend
simply translates to instantiate the PolarForm
structure on the GPU: polar = PolarForm("case9.m", CUDABackend())
.
Interface
All functions are following AutoDiff.AbstractExpression
's interface. The structure of the network is specified by the PolarForm
we pass as an argument in the constructor. For instance, we build a new PolarBasis
expression associated to case9
directly as
julia> basis = ExaPF.PolarBasis(polar)
PolarBasis (AbstractExpression)
Each expression as a given dimension, given by
julia> length(basis)
27
In ExaPF, the inputs and the parameters are stored inside a NetworkStack
structure. Evaluating the basis $\psi$ naturally translates to
julia> basis(stack)
27-element Vector{Float64}:
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
0.0
⋮
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
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
:
julia> psi = zeros(length(basis)) ;
julia> basis(psi, stack);
Compose expressions together
In the LTMR2020 model, the polar basis $\psi(v, \theta)$ depends only on the voltage magnitudes and the voltage angles. However, this is not the case for the other functions (power flow balance, line flows, ...), which all depends on the basis $\psi(v, \theta)$.
In ExaPF, one has to build manually the vectorized expression tree associated to the power flow model. Luckily, evaluating the LTMR2020 simply amounts to compose functions together with the polar basis $\psi(v, \theta)$. ExaPF overloads the function ∘
to compose functions with a PolarBasis
instance. The power flow balance can be evaluated as
julia> pflow = ExaPF.PowerFlowBalance(polar) ∘ basis;
which returns a ComposedExpressions
structure.
The function pflow
follows the same API, as any regular AutoDiff.AbstractExpression
.
julia> n_balance = length(pflow)
14
julia> round.(pflow(stack); digits=6) # evaluate the power flow balance
14-element Vector{Float64}:
-1.63
-0.85
0.0
0.9
0.0
1.0
0.0
1.25
-0.167
0.042
-0.2835
0.171
-0.2275
0.259
When we evaluate a ComposedExpressions
, ExaPF first computes the basis $\psi(v, \theta)$ inside stack.psi
, and then ExaPF uses the values in stack.psi
to evaluate the final result.
The procedure remains the same if one wants to evaluate the LineFlows
or the PowerGenerationBounds
. For instance, evaluating the line flows amounts to
julia> line_flows = ExaPF.LineFlows(polar) ∘ basis;
julia> line_flows(stack)
18-element Vector{Float64}:
0.0
0.006241000000000099
0.0320410000000001
0.0
0.010920249999999961
0.005550250000000068
0.0
0.02340899999999987
0.007743999999999858
0.0
0.006241000000000099
0.0320410000000001
0.0
0.010920249999999961
0.005550250000000068
0.0
0.02340899999999987
0.007743999999999858