1 - ⚡ Zap Basics ⚡¶
In this notebook, we introduce the basic types and functionality in Zap. We will manually construct a small electricity network, solve a dispatch problem using both CVXPY and ADMM, and analyze the results.
Creating a Network
Solving Dispatch Problems
Analyzing Results
[1]:
import zap
import torch
import numpy as np
import cvxpy as cp
Creating a Network¶
Networks in Zap consist of a Network object, which defines the underlying network, and a list of Device objects, which are electrical components attached to the network. Devices include generators, loads, transmission lines, batteries, and all other components that connect to the network. Devices can also be used to model financial contracts, such as bids and offers for energy or transmission congestion contracts.
Constructing Devices¶
First, let’s initialize a simple network with 3 nodes.
[2]:
net = zap.PowerNetwork(num_nodes=3)
Now, let’s add a load to the network.
[3]:
load = zap.Load(
num_nodes=net.num_nodes,
# The terminal of the device is the node to which it connects
terminal=np.array([0]),
# The load argument is a time series of the power consumed by the device
# This is 2-d array, see more below!
load=np.array([[10.0, 15.0, 20.0]]),
# This is the curtailment cost of the device, which we will assume is $500 / MWh
linear_cost=np.array([500.0]),
)
All data passed to the device must be numpy arrays or of a similar type (e.g., a PyTorch Tensor). We can also create groups of devices at the same time. Let’s add two generators to the network.
[4]:
generators = zap.Generator(
num_nodes=net.num_nodes,
# Since we have two generators, we specify two terminals
terminal=np.array([1, 2]),
# Nominal capacity refers to the generator nameplate capacity
nominal_capacity=np.array([50.0, 25.0]),
# This is the marginal cost of generation
linear_cost=np.array([0.1, 30.0]),
# Emissions rates are optional, but useful for many applications
emission_rates=np.array([0.0, 500.0]),
# Dynamic capacity refers to the time-varying capacity factor of each generator
dynamic_capacity=np.array(
[
[0.1, 0.5, 0.1],
[1.0, 1.0, 1.0],
]
),
)
Since we have two generators, we specify both their parameters at the same time. Some of the parameters, such as costs or dynamic capacities, can also time varying, in which case the user must pass a 2d-array of data. The first dimension of the data specifies which device is being referred to, and the second dimension of the data specifies the time period. If you pass a 1d-array, Zap will automatically assume this is a static quantity that does not change over time.
Finally, let’s initialize a few transmission lines to connect our network.
[5]:
lines = zap.ACLine(
num_nodes=net.num_nodes,
# AC lines are two-terminal devices, so we specify both a source and sink terminal for each line
# We will build the classic 3-bus "triangle" network
source_terminal=np.array([0, 1, 2]),
sink_terminal=np.array([1, 2, 0]),
nominal_capacity=np.array([10.0, 20.0, 30.0]),
# This is the per-MW susceptance, so the susceptance of each line is its nominal capacity
# times its susceptance
# We will give the lines uniform total susceptance
susceptance=np.array([1 / 10.0, 1 / 20.0, 1 / 30.0]),
# Lines can also have time-vary capacities to simulate outages or de-ratings
# In this example, we will just make them static
capacity=np.ones(3),
)
We will also add an electrical ground device to specify which node is the reference bus.
[6]:
ground = zap.Ground(num_nodes=net.num_nodes, terminal=np.array([0]))
Putting the Network Together¶
To finish up, we will create a list of device (groups). Together with the network and the time horizon, this list fully specifies an electrical network!
[7]:
devices = [load, generators, lines, ground]
time_horizon = 3
net, devices, time_horizon # Fully network specification
[7]:
(PowerNetwork(num_nodes=3),
[Load(num_nodes=3, terminal=array([0]), nominal_capacity=array([[1.]]), min_power=array([[-10., -15., -20.]]), max_power=array([[0., 0., 0.]]), capital_cost=None, emission_rates=None, load=array([[10., 15., 20.]]), linear_cost=array([[500.]]), quadratic_cost=None),
Generator(num_nodes=3, terminal=array([1, 2]), nominal_capacity=array([[50.],
[25.]]), min_power=array([[0., 0., 0.],
[0., 0., 0.]]), max_power=array([[0.1, 0.5, 0.1],
[1. , 1. , 1. ]]), dynamic_capacity=array([[0.1, 0.5, 0.1],
[1. , 1. , 1. ]]), linear_cost=array([[ 0.1],
[30. ]]), quadratic_cost=None, capital_cost=None, emission_rates=array([[ 0.],
[500.]]), min_nominal_capacity=None, max_nominal_capacity=None),
ACLine(num_nodes=3, source_terminal=array([0, 1, 2]), sink_terminal=array([1, 2, 0]), min_power=array([[-1.],
[-1.],
[-1.]]), max_power=array([[1.],
[1.],
[1.]]), linear_cost=array([[0.],
[0.],
[0.]]), quadratic_cost=None, nominal_capacity=array([[10.],
[20.],
[30.]]), capital_cost=None, slack=0.0, min_nominal_capacity=None, max_nominal_capacity=None, reconductoring_cost=None, reconductoring_threshold=None),
Ground(num_nodes=3, terminal=array([0]), voltage=array([[0.]]))],
3)
Before we move on, here are a couple of notes about constructing devices.
Many arguments (capital cost, emissions rates, line slacks, etc.) are optional and introduce additional functionality, e.g., for planning studies.
For performance reasons, it’s best to specify a few device groups (ideally, one for each type of device) with many devices per group, instead of many groups with one or just a few devices per group. This is so that computations can be efficiently vectorized, both on CPU and GPU machines.
Remember that the first dimension (rows) of a parameter specifies which device in the group the parameter is for, and the second dimension (columns) specifies the time period. Dynamic data is always 2-dimensional!
Solving Dispatch Problems¶
Now that we’ve defined an electricity dispatch problem, let’s solve using two different methods.
cvxpy- This method will build a model incvxpyand send it to an off-the-shelf solver, such as ECOS or Clarabel. You can also use commerical solvers like Mosek. This approach is best for small to medium problems and finds highly accurate solutions.admm- This method will use a PyTorch implementation of the alternating direction method of multipliers (ADMM). This can be run on either a CPU or GPU and scales well to larger problems. In general, this method is only capable of finding medium accuracy solutions within a reasonable amount of time.
Solving with CVXPY¶
To solve with cvxpy, we simply call net.dispatch on our devices and (optionally) specify a solver.
[8]:
outcome = net.dispatch(devices, time_horizon, solver=cp.CLARABEL, add_ground=False)
outcome
[8]:
DispatchOutcome(phase_duals=[None, None, [array([[-7.21619994e-08, -1.48353913e+01, -5.69716112e-07],
[-7.21620029e-08, -1.48353913e+01, -5.69716118e-07],
[-7.21619968e-08, -1.48353913e+01, -5.69716103e-07]]), array([[7.21620016e-08, 1.48353913e+01, 5.69716116e-07],
[7.21619999e-08, 1.48353913e+01, 5.69716109e-07],
[7.21619944e-08, 1.48353913e+01, 5.69716097e-07]])], [array([[ 5.02100773e-15, -2.31869585e-15, 1.46905872e-14]])]], local_equality_duals=[[], [], [array([[29.99999974, 29.77078293, 30.00000068],
[29.99999942, 0.10000012, 29.99999949],
[29.9999996 , 14.93539156, 30.00000011]]), array([[7.21619994e-08, 1.48353913e+01, 5.69716112e-07],
[7.21620041e-08, 1.48353913e+01, 5.69716121e-07],
[7.21619944e-08, 1.48353913e+01, 5.69716097e-07]])], [array([[29.99999974, 29.77078293, 30.00000068]]), array([[ 5.02100773e-15, -2.31869585e-15, 1.46905872e-14]])]], local_inequality_duals=[[array([[470.0000004 , 470.22921717, 469.99999943]]), array([[1.40746090e-07, 9.60038240e-08, 1.09352091e-07]])], [array([[4.28799138e-07, 1.78300241e-07, 3.37175404e-07],
[4.56793016e-07, 1.50646085e+01, 1.61228178e-07]]), array([[2.98999998e+01, 3.03066738e-07, 2.98999998e+01],
[6.01385815e-08, 4.00736531e-08, 2.74690419e-07]])], [array([[4.04079311e-07, 4.45061741e+01, 1.83415172e-06],
[7.31746157e-08, 6.55602650e-08, 9.83955076e-08],
[5.59249063e-08, 5.71740987e-08, 6.12132993e-08]]), array([[9.44655495e-08, 7.82423597e-08, 9.86035365e-08],
[1.17743892e-07, 1.51057653e-07, 1.05392772e-07],
[9.05760184e-08, 7.98719289e-08, 6.55684197e-08]])], []], local_variables=[None, None, None, None], power=[[array([[-10., -15., -20.]])], [array([[5.00000000e+00, 1.50000000e+01, 4.99999999e+00],
[5.00000001e+00, 9.03717974e-09, 1.50000000e+01]])], [array([[ 5.00000000e+00, 1.00000000e+01, 8.33333333e+00],
[ 3.64665045e-09, -4.99999999e+00, 3.33333334e+00],
[-5.00000000e+00, -5.00000000e+00, -1.16666667e+01]]), array([[-5.00000000e+00, -1.00000000e+01, -8.33333333e+00],
[-3.64607602e-09, 4.99999999e+00, -3.33333334e+00],
[ 5.00000000e+00, 5.00000000e+00, 1.16666667e+01]])], [array([[ 2.59944050e-13, -1.01386661e-13, -6.75087621e-13]])]], angle=[None, None, [array([[-7.14401061e-14, 8.38978856e-14, -5.64017010e-13],
[ 5.00000000e+00, 1.00000000e+01, 8.33333333e+00],
[ 5.00000000e+00, 5.00000000e+00, 1.16666667e+01]]), array([[ 5.00000000e+00, 1.00000000e+01, 8.33333333e+00],
[ 5.00000000e+00, 5.00000000e+00, 1.16666667e+01],
[ 7.14401067e-14, -8.38978841e-14, 5.64017011e-13]])], [array([[1.19208492e-22, 3.02650078e-22, 2.33199737e-22]])]], prices=array([[29.99999974, 29.77078293, 30.00000068],
[29.99999942, 0.10000012, 29.99999949],
[29.9999996 , 14.93539156, 30.00000011]]), global_angle=array([[2.38416984e-22, 6.05300156e-22, 4.66399474e-22],
[5.00000000e+00, 1.00000000e+01, 8.33333333e+00],
[5.00000000e+00, 5.00000000e+00, 1.16666667e+01]]), problem=Problem(Minimize(Expression(AFFINE, UNKNOWN, ())), [Equality(Expression(AFFINE, UNKNOWN, (3, 3)), Constant(CONSTANT, ZERO, ())), Equality(Expression(AFFINE, UNKNOWN, (3, 3)), Variable((3, 3), var7)), Equality(Expression(AFFINE, UNKNOWN, (3, 3)), Variable((3, 3), var8)), Equality(Expression(AFFINE, UNKNOWN, (1, 3)), Variable((1, 3), var9)), Equality(Expression(AFFINE, UNKNOWN, (3, 3)), Constant(CONSTANT, ZERO, ())), Equality(Expression(AFFINE, UNKNOWN, (3, 3)), Constant(CONSTANT, ZERO, ())), Equality(Variable((1, 3), var6), Constant(CONSTANT, ZERO, ())), Equality(Expression(AFFINE, UNKNOWN, (1, 3)), Constant(CONSTANT, ZERO, ())), Inequality(Expression(AFFINE, UNKNOWN, (1, 3))), Inequality(Expression(AFFINE, UNKNOWN, (1, 3))), Inequality(Expression(AFFINE, UNKNOWN, (2, 3))), Inequality(Expression(AFFINE, UNKNOWN, (2, 3))), Inequality(Expression(AFFINE, UNKNOWN, (3, 3))), Inequality(Expression(AFFINE, UNKNOWN, (3, 3)))]), ground=None)
As you can see, the result of a grid dispatch is a complicated object. We will address that in a second.
Solving with ADMM¶
Solving with ADMM is a little more complicated and requires two steps:
Transfering device data to PyTorch
Initializing an ADMM solver object.
[9]:
machine = "cpu" # Pick "cuda" for a machine with an Nvidia GPU
dtype = torch.float32
[10]:
from zap.admm import ADMMSolver
[11]:
admm_devices = [d.torchify(machine=machine, dtype=dtype) for d in devices]
[12]:
admm = ADMMSolver(
machine=machine,
dtype=dtype,
atol=1e-6,
rtol=1e-6,
)
solution_admm, history_admm = admm.solve(net, admm_devices, time_horizon)
ADMM converged in 468 iterations.
[13]:
solution_admm.power
[13]:
[[tensor([[-10., -15., -20.]])],
[tensor([[ 5.0000, 15.0000, 5.0000],
[ 5.0001, 0.0000, 15.0001]])],
[tensor([[ 5.0000e+00, 1.0000e+01, 8.3333e+00],
[ 2.9945e-05, -5.0000e+00, 3.3334e+00],
[-5.0000e+00, -5.0000e+00, -1.1667e+01]]),
tensor([[-5.0000e+00, -1.0000e+01, -8.3333e+00],
[-2.9945e-05, 5.0000e+00, -3.3334e+00],
[ 5.0000e+00, 5.0000e+00, 1.1667e+01]])],
[tensor([[0., 0., 0.]])]]
[14]:
# ADMM solutions need to be cast to a standard DispatchOutcome
outcome_admm = solution_admm.as_outcome()
outcome_admm
[14]:
DispatchOutcome(phase_duals=[None, None, [tensor([[-3.3689e-05, 3.5381e+00, 3.5673e-05],
[-2.6349e-05, 3.5381e+00, 5.4371e-05],
[ 1.9117e-05, 3.5381e+00, -2.1447e-05]]), tensor([[ 1.8560e-05, -3.5381e+00, -1.4193e-05],
[-1.1102e-05, -3.5381e+00, 9.2870e-06],
[-4.8603e-05, -3.5381e+00, 7.4019e-05]])], [tensor([[ 8.0933e-05, 2.5490e-05, -1.1069e-04]])]], local_equality_duals=None, local_inequality_duals=None, local_variables=[None, None, None, None], power=[[tensor([[-10., -15., -20.]])], [tensor([[ 5.0000, 15.0000, 5.0000],
[ 5.0001, 0.0000, 15.0001]])], [tensor([[ 5.0000e+00, 1.0000e+01, 8.3333e+00],
[ 2.9945e-05, -5.0000e+00, 3.3334e+00],
[-5.0000e+00, -5.0000e+00, -1.1667e+01]]), tensor([[-5.0000e+00, -1.0000e+01, -8.3333e+00],
[-2.9945e-05, 5.0000e+00, -3.3334e+00],
[ 5.0000e+00, 5.0000e+00, 1.1667e+01]])], [tensor([[0., 0., 0.]])]], angle=[None, None, [tensor([[1.1444e-05, 0.0000e+00, 1.2398e-05],
[5.0000e+00, 1.0000e+01, 8.3334e+00],
[5.0001e+00, 5.0000e+00, 1.1667e+01]]), tensor([[5.0000e+00, 1.0000e+01, 8.3334e+00],
[5.0001e+00, 5.0000e+00, 1.1667e+01],
[1.2398e-05, 0.0000e+00, 1.3351e-05]])], [tensor([[0., 0., 0.]])]], prices=tensor([[30.0001, 7.1762, 29.9998],
[30.0001, 0.1000, 29.9998],
[30.0000, 3.6381, 29.9999]]), global_angle=None, problem=None, ground=None)
Analyzing Results¶
Results are packaged into a hierachically structured DispatchOutcome object.
At the top level, a
DispatchOutcomehas several fields:power,angle,prices, and a few others. You can access these fields like any other Python field, e.g.,outcome.power.Each field contains either device-specific information or global information. Device-specific fields will contain a list of length
len(devices), where theith entry in the list contains information specific to theith device. Global fields contain a 2d array of size(num_nodes, time_horizon). You can access the information for deviceiusing normal indexing, e.g.,outcome.power[i].For device-specific information, each block of information is further broken down by the terminal of the device. Many devices, such as generator, loads, and batteries, have just a single terminal and will always be indexed as
outcome.power[i][0]. Some devices, like transmission lines, have two or more terminals. The data for terminaljis stored inoutcome.power[i][j].Finally, the data for terminal
jof deviceiis just a 2d array of size(num_devices, time_horizon), wherenum_devices = devices[i].num_devicesis the number of devices in device groupi.
[15]:
outcome.prices # Global
[15]:
array([[29.99999974, 29.77078293, 30.00000068],
[29.99999942, 0.10000012, 29.99999949],
[29.9999996 , 14.93539156, 30.00000011]])
[16]:
outcome.power # Device-specific
[16]:
[[array([[-10., -15., -20.]])],
[array([[5.00000000e+00, 1.50000000e+01, 4.99999999e+00],
[5.00000001e+00, 9.03717974e-09, 1.50000000e+01]])],
[array([[ 5.00000000e+00, 1.00000000e+01, 8.33333333e+00],
[ 3.64665045e-09, -4.99999999e+00, 3.33333334e+00],
[-5.00000000e+00, -5.00000000e+00, -1.16666667e+01]]),
array([[-5.00000000e+00, -1.00000000e+01, -8.33333333e+00],
[-3.64607602e-09, 4.99999999e+00, -3.33333334e+00],
[ 5.00000000e+00, 5.00000000e+00, 1.16666667e+01]])],
[array([[ 2.59944050e-13, -1.01386661e-13, -6.75087621e-13]])]]
[17]:
outcome.power[2] # Power flows for device 2, which we made the transmission lines
[17]:
[array([[ 5.00000000e+00, 1.00000000e+01, 8.33333333e+00],
[ 3.64665045e-09, -4.99999999e+00, 3.33333334e+00],
[-5.00000000e+00, -5.00000000e+00, -1.16666667e+01]]),
array([[-5.00000000e+00, -1.00000000e+01, -8.33333333e+00],
[-3.64607602e-09, 4.99999999e+00, -3.33333334e+00],
[ 5.00000000e+00, 5.00000000e+00, 1.16666667e+01]])]
[18]:
# Indexing the second terminal of the lines---the result is an object of size (num_ac_lines, time_horizon)
outcome.power[2][1]
[18]:
array([[-5.00000000e+00, -1.00000000e+01, -8.33333333e+00],
[-3.64607602e-09, 4.99999999e+00, -3.33333334e+00],
[ 5.00000000e+00, 5.00000000e+00, 1.16666667e+01]])
[19]:
# Indexing the first (and only) terminal of the generators
outcome.power[1][0]
[19]:
array([[5.00000000e+00, 1.50000000e+01, 4.99999999e+00],
[5.00000001e+00, 9.03717974e-09, 1.50000000e+01]])
[20]:
# Power output of generator 0 at timestep 2
outcome.power[1][0][0, 2]
[20]:
4.999999989280596