AbstractVariable

Variables corresponding to a particular formulation.

AbstractFormulation

Interface between the data and the mathemical formulation.

State <: AbstractVariable

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

Control <: AbstractVariable

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

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 backend (CPU() and CUDABackend() are currently supported).

Example

julia> const PS = ExaPF.PowerSystem;

julia> network_data = PS.load_case("case9.m");

julia> polar = PolarForm(network_data, ExaPF.CPU())
Polar formulation (instantiated on backend CPU(false))
Network characteristics:
    #buses:      9  (#slack: 1  #PV: 2  #PQ: 6)
    #generators: 3
    #lines:      9
giving a mathematical formulation with:
    #controls:   5
    #states  :   14

BlockPolarForm{T, IT, VT, MT} <: AbstractFormulation

Block polar formulation: duplicates k different polar models to evaluate them in parallel.

load_polar(case, backend=CPU(); dir=PS.EXADATA)

Load a PolarForm instance from the specified benchmark library dir on the target backend (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 backend CPU(false))
Network characteristics:
    #buses:      9  (#slack: 1  #PV: 2  #PQ: 6)
    #generators: 3
    #lines:      9
giving a mathematical formulation with:
    #controls:   5
    #states  :   14

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

Examples

julia> polar = ExaPF.load_polar("case9");

julia> stack = ExaPF.NetworkStack(polar)
21-elements NetworkStack{Vector{Float64}}

julia> stack.vmag
9-element Vector{Float64}:
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0

init!(polar::PolarForm, stack::NetworkStack)

Set stack.input with the initial values specified in the base PS.PowerNetwork object.

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}:
  1
  2
  3
 20
 21

mapping(polar::PolarForm, ::State)

Return the mapping associated to the State() in NetworkStack according to the polar formulation PolarForm.

Examples

julia> polar = ExaPF.load_polar("case9");

julia> mapu = ExaPF.mapping(polar, State())
14-element Vector{Int64}:
 11
 12
 13
 14
 15
 16
 17
 18
  4
  5
  6
  7
  8
  9

PowerFlowProblem

A mutable struct representing a power flow problem.

Fields

• form::AbstractFormulation: The problem formulation (e.g., PolarForm, BlockPolarForm)
• stack::AbstractNetworkStack: Stack containing network variables and parameters
• powerflow::AD.AbstractExpression: Power flow balance expression
• linear_solver::LS.AbstractLinearSolver: Linear solver for Newton-Raphson iterations
• non_linear_solver::AbstractNonLinearSolver: Non-linear solver configuration
• mapx::Vector{Int}: Mapping indices for state variables
• jac::AD.AbstractJacobian: Jacobian matrix for automatic differentiation
• conv::ConvergenceStatus: Convergence status of the solver
• backend::KA.Backend: Computation backend (CPU or GPU)

See also

run_pf, solve!

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");

julia> stack = ExaPF.NetworkStack(polar);

julia> conv = run_pf(polar, stack);

julia> conv.has_converged
true

julia> conv.n_iterations
4

solve!(prob::PowerFlowProblem)

Re-solve an existing power flow problem.

This function re-runs the non-linear solver on the problem, which is useful after modifying problem parameters (e.g., via set_active_load! or set_reactive_load!).

Arguments

• prob::PowerFlowProblem: The power flow problem to solve

Returns

• Convergence status of the solver

Examples

# Create and solve initial problem
prob = PowerFlowProblem("case9.m", CPU(), :polar)
solve!(prob)

# Modify parameters and re-solve
set_active_load!(prob, ones(9) * 0.8)
solve!(prob)

See also

run_pf, set_active_load!, set_reactive_load!

nlsolve!(
    algo::NewtonRaphson,
    jac::Jacobian,
    stack::NetworkStack;
    linear_solver=DirectSolver(jac.J),
    nl_buffer=NLBuffer(size(jac, 2)),
)

Solve the nonlinear system of equations $g(x) = 0$ with a NewtonRaphson algorithm. At each iteration, we update the variable $x$ as

\[ x_{k+1} = x_{k} - (∇g_k)^{-1} g(x_k) \]

till $\| g(x_k) \| < ε_{tol}$

In the implementation,

• the function $g$ is specified in jac.func,
• the initial variable $x_0$ in stack::NetworkStack (with mapping jac.map),
• the Jacobian $∇g$ is computed automatically in jac, with automatic differentiation.

Note that stack is modified inplace during the iterations of algorithm.

The Jacobian jac should be instantied before calling this function. By default, the linear system $(∇g_k)^{-1} g(x_k)$ is solved using a LU factorization. You can specify a different linear solver by changing the optional argument linear_solver.

Arguments

• algo::NewtonRaphon: Newton-Raphson object, storing the options of the algorithm
• jac::Jacobian: Stores the function $g$ and its Jacobian $∇g$. The Jacobian is updated with automatic differentiation.
• stack::NetworkStack: initial values
• linear_solver::AbstractLinearSolver: linear solver used to compute the Newton step
• nl_buffer::NLBuffer: buffer storing the residual vector and the descent direction Δx. Can be reused to avoid unecessary allocations.

Examples

julia> polar = ExaPF.load_polar("case9");

julia> powerflow = ExaPF.PowerFlowBalance(polar) ∘ ExaPF.Basis(polar);

julia> jx = ExaPF.Jacobian(polar, powerflow, State());

julia> stack = ExaPF.NetworkStack(polar);

julia> conv = ExaPF.nlsolve!(NewtonRaphson(), jx, stack);

julia> conv.has_converged
true

julia> conv.n_iterations
4

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

get_active_load(prob::PowerFlowProblem)

Get the active power demand (Pd) values from the power flow problem.

Arguments

• prob::PowerFlowProblem: The power flow problem instance

Returns

• Vector: Active power demand values for all buses in the system

See also

get_reactive_load, set_active_load!

get_reactive_load(prob::PowerFlowProblem)

Get the reactive power demand (Qd) values from the power flow problem.

Arguments

• prob::PowerFlowProblem: The power flow problem instance

Returns

• Vector: Reactive power demand values for all buses in the system

See also

get_active_load, set_reactive_load!

set_active_load!(prob::PowerFlowProblem, pd::Vector{Float64})

Set the active power demand (Pd) values for the power flow problem.

This function modifies the problem in-place, updating the active power demand parameters in the network stack.

Arguments

• prob::PowerFlowProblem: The power flow problem instance to modify
• pd::Vector{Float64}: New active power demand values for all buses

Returns

• PowerFlowProblem: The modified problem instance

Throws

• AssertionError: If the length of pd does not match the number of buses

Examples

prob = PowerFlowProblem("case9.m", CPU(), :polar)
pd_new = ones(9) * 0.5  # Set all buses to 0.5 p.u.
set_active_load!(prob, pd_new)

See also

set_reactive_load!, get_active_load

set_reactive_load!(prob::PowerFlowProblem, qd::Vector{Float64})

Set the reactive power demand (Qd) values for the power flow problem.

This function modifies the problem in-place, updating the reactive power demand parameters in the network stack.

Arguments

• prob::PowerFlowProblem: The power flow problem instance to modify
• qd::Vector{Float64}: New reactive power demand values for all buses

Returns

• PowerFlowProblem: The modified problem instance

Throws

• AssertionError: If the length of qd does not match the number of buses

Examples

prob = PowerFlowProblem("case9.m", CPU(), :polar)
qd_new = ones(9) * 0.2  # Set all buses to 0.2 p.u.
set_reactive_load!(prob, qd_new)

See also

set_active_load!, get_reactive_load

get_voltage_magnitude(prob::PowerFlowProblem)

Get the voltage magnitude values from the power flow problem solution.

Arguments

• prob::PowerFlowProblem: The power flow problem instance

Returns

• Vector: Voltage magnitude values (in per-unit) for all buses in the system

See also

get_voltage_angle, get_solution

get_voltage_angle(prob::PowerFlowProblem)

Get the voltage angle values from the power flow problem solution.

Arguments

• prob::PowerFlowProblem: The power flow problem instance

Returns

• Vector: Voltage angle values (in radians) for all buses in the system

See also

get_voltage_magnitude, get_solution

get_solution(prob::PowerFlowProblem)

Get the complete solution vector from the power flow problem.

This function returns the state variables (voltage angles and magnitudes) that were solved by the Newton-Raphson method, in the order specified by the mapping indices.

Arguments

• prob::PowerFlowProblem: The power flow problem instance

Returns

• Vector: Solution vector containing the state variables

See also

get_voltage_angle, get_voltage_magnitude, solve!

Basis{VI, MT} <: AbstractExpression
Basis(polar::AbstractPolarFormulation)

Implement the LKMR nonlinear basis. Takes as input the voltage magnitudes vmag and the voltage angles vang and returns

\[ \begin{aligned} & \psi_\ell^C(v, \theta) = v^f v^t \cos(\theta_f - \theta_t) \quad \forall \ell = 1, \cdots, n_\ell \\ & \psi_\ell^S(v, \theta) = v^f v^t \sin(\theta_f - \theta_t) \quad \forall \ell = 1, \cdots, n_\ell \\ & \psi_k(v, \theta) = v_k^2 \quad \forall k = 1, \cdots, n_b \end{aligned}\]

Dimension: 2 * n_lines + n_bus

Complexity

3 n_lines + n_bus mul, n_lines cos and n_lines sin

Examples

julia> polar = ExaPF.load_polar("case9");

julia> stack = ExaPF.NetworkStack(polar);

julia> basis = ExaPF.Basis(polar)
Basis (AbstractExpression)

julia> basis(stack)
27-element Vector{Float64}:
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 0.0
 ⋮
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0

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\} \\ q_i &= v_i \sum_{j}^{n} v_j (g_{ij}\sin{(\theta_i - \theta_j)} - b_{ij}\cos{(\theta_i - \theta_j})) \,. & ∀ i ∈ \{PQ\} \end{aligned}\]

Require composition with Basis.

Dimension: n_pv + 2 * n_pq

Complexity

2 SpMV

Examples

julia> polar = ExaPF.load_polar("case9");

julia> stack = ExaPF.NetworkStack(polar);

julia> powerflow = ExaPF.PowerFlowBalance(polar) ∘ ExaPF.Basis(polar);

julia> round.(powerflow(stack); digits=6)
14-element Vector{Float64}:
 -1.63
 -0.85
  0.0
  0.9
  0.0
  1.0
  0.0
  1.25
 -0.167
  0.042
 -0.2835
  0.171
 -0.2275
  0.259

julia> run_pf(polar, stack); # solve powerflow equations

julia> isapprox(powerflow(stack), zeros(14); atol=1e-8)
true

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);

julia> voltage_pq = ExaPF.VoltageMagnitudeBounds(polar);

julia> voltage_pq(stack)
6-element Vector{Float64}:
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0

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

Dimension: n_pv + 2 n_ref

Complexity

1 copyto, 1 SpMV

Examples

julia> polar = ExaPF.load_polar("case9");

julia> stack = ExaPF.NetworkStack(polar);

julia> run_pf(polar, stack); # solve powerflow equations

julia> power_generators = ExaPF.PowerGenerationBounds(polar) ∘ ExaPF.Basis(polar);

julia> round.(power_generators(stack); digits=6)
4-element Vector{Float64}:
  0.719547
  0.24069
  0.144601
 -0.03649

LineFlows{VT, MT}
LineFlows(polar)

Implement thermal limit constraints on the lines of the network.

Require composition with Basis.

Dimension: 2 * n_lines

Complexity

4 SpMV, 4 * n_lines quadratic, 2 * n_lines add

Examples

julia> polar = ExaPF.load_polar("case9");

julia> stack = ExaPF.NetworkStack(polar);

julia> run_pf(polar, stack); # solve powerflow equations

julia> line_flows = ExaPF.LineFlows(polar) ∘ ExaPF.Basis(polar);

julia> round.(line_flows(stack); digits=6)
18-element Vector{Float64}:
 0.575679
 0.094457
 0.379983
 0.723832
 0.060169
 0.588673
 2.657418
 0.748943
 0.295351
 0.560817
 0.112095
 0.38625
 0.728726
 0.117191
 0.585164
 2.67781
 0.726668
 0.215497

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 Basis to evaluate the cost of the reference generator.

Dimension: 1

Complexity

1 SpMV, 1 sum

Examples

julia> polar = ExaPF.load_polar("case9");

julia> stack = ExaPF.NetworkStack(polar);

julia> cost = ExaPF.CostFunction(polar) ∘ ExaPF.Basis(polar);

julia> cost(stack)
1-element Vector{Float64}:
 4509.0275

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)]

ComposedExpressions{Expr1<:Basis, 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 Basis is supported for expr1.