pandapower now has an interface to PowerModels.jl that can be used for efficient power system optimization.
Because pandapower:
So using pandapower to define the grid models and then using PowerModels for the optimization really gives you the best of all worlds - you can use the rich environment of Python libraries, the sophisticated element models of pandapower, the modular optimization framework of PowerModels and the efficient mathematical modeling of JuMP.
So here is an example of how it works. First, we create a grid in pandapower. Here, we create a meshed 110kV grid with four buses that is fed from an 220kV network through a 3-Winding transformer.
import pandapower as pp
import numpy as np
net = pp.create_empty_network()
min_vm_pu = 0.95
max_vm_pu = 1.05
#create buses
bus1 = pp.create_bus(net, vn_kv=220., geodata=(5,9), min_vm_pu=min_vm_pu, max_vm_pu=max_vm_pu)
bus2 = pp.create_bus(net, vn_kv=110., geodata=(6,10), min_vm_pu=min_vm_pu, max_vm_pu=max_vm_pu)
bus3 = pp.create_bus(net, vn_kv=110., geodata=(10,9), min_vm_pu=min_vm_pu, max_vm_pu=max_vm_pu)
bus4 = pp.create_bus(net, vn_kv=110., geodata=(8,8), min_vm_pu=min_vm_pu, max_vm_pu=max_vm_pu)
bus5 = pp.create_bus(net, vn_kv=110., geodata=(6,8), min_vm_pu=min_vm_pu, max_vm_pu=max_vm_pu)
#create 220/110/110 kV 3W-transformer
pp.create_transformer3w_from_parameters(net, bus1, bus2, bus5, vn_hv_kv=220, vn_mv_kv=110,
vn_lv_kv=110, vk_hv_percent=10., vk_mv_percent=10.,
vk_lv_percent=10., vkr_hv_percent=0.5,
vkr_mv_percent=0.5, vkr_lv_percent=0.5, pfe_kw=10,
i0_percent=0.1, shift_mv_degree=0, shift_lv_degree=0,
sn_hv_mva=100, sn_mv_mva=50, sn_lv_mva=50)
#create 110 kV lines
l1 = pp.create_line(net, bus2, bus3, length_km=70., std_type='149-AL1/24-ST1A 110.0')
l2 = pp.create_line(net, bus3, bus4, length_km=50., std_type='149-AL1/24-ST1A 110.0')
l3 = pp.create_line(net, bus4, bus2, length_km=40., std_type='149-AL1/24-ST1A 110.0')
l4 = pp.create_line(net, bus4, bus5, length_km=30., std_type='149-AL1/24-ST1A 110.0')
#create loads
pp.create_load(net, bus2, p_mw=60)
pp.create_load(net, bus3, p_mw=70)
pp.create_load(net, bus4, p_mw=10)
#create generators
g1 = pp.create_gen(net, bus1, p_mw=40, min_p_mw=0, max_p_mw=200, vm_pu=1.01, slack=True)
pp.create_poly_cost(net, g1, 'gen', cp1_eur_per_mw=1)
g2 = pp.create_gen(net, bus3, p_mw=40, min_p_mw=0, max_p_mw=200, vm_pu=1.01)
pp.create_poly_cost(net, g2, 'gen', cp1_eur_per_mw=3)
g3 = pp.create_gen(net, bus4, p_mw=50, min_p_mw=0, max_p_mw=200, vm_pu=1.01)
pp.create_poly_cost(net, g3, 'gen', cp1_eur_per_mw=3)
net
This pandapower network includes the following parameter tables: - bus (5 elements) - load (3 elements) - gen (3 elements) - line (4 elements) - trafo3w (1 element) - poly_cost (3 elements) - bus_geodata (5 elements)
Note that PowerModels does not have a 3W-transformer model, but since pandapower includes the equations to calculates the equivalent branches for the 3W-transformers, it is possible to optimize grids with 3W-transformers in PowerModels through the pandapower interface. The same is true for other complex transformer models, switches/breaker, extended ward equivalents etc.
Let's have a look at the grid we created with pandapowers plotting module:
import pandapower.plotting as plot
%matplotlib inline
plot.simple_plot(net)
<matplotlib.axes._subplots.AxesSubplot at 0xaebf630>
Now lets run an OPF through PowerModels and look at the results (Note that the first time the runpm function is called, Julia is started in the background, which may take some time):
pp.runpm_ac_opf(net)
Since Generator 1 has the lowest cost, all required power is supplied through this generator:
net.res_gen
p_mw | q_mvar | va_degree | vm_pu | |
---|---|---|---|---|
0 | 144.851017 | 11.037217 | 7.117490e-27 | 1.050000 |
1 | 0.000000 | 6.055925 | -1.680643e+01 | 1.013179 |
2 | 0.000000 | 13.126070 | -1.298596e+01 | 1.040498 |
This however leeds to an overload in the three-winding transformer, through which g1 is connected:
net.res_trafo3w.loading_percent
0 190.835372 Name: loading_percent, dtype: float64
Let's set some constraints for the 3W-transformer and the lines and rerun the OPF:
net.trafo3w["max_loading_percent"] = 50
net.line["max_loading_percent"] = 20
pp.runpm_ac_opf(net)
The constraints are complied with for all lines and the 3W transformer:
net.res_trafo3w.loading_percent
0 49.97477 Name: loading_percent, dtype: float64
net.res_line.loading_percent
0 19.047619 1 12.295418 2 19.207577 3 7.088595 Name: loading_percent, dtype: float64
The power is now generated by a mixture of the generators:
net.res_gen
p_mw | q_mvar | va_degree | vm_pu | |
---|---|---|---|---|
0 | 24.531626 | -6.267018 | 2.961294e-28 | 1.031058 |
1 | 98.101643 | -7.928705 | -3.576364e-01 | 1.050000 |
2 | 18.054898 | 9.928328 | -1.446870e+00 | 1.048757 |
Apart from the AC OPF used in the example above, pandapower also has an interface to run the DC OPF:
pp.runpm_dc_opf(net)
net.res_bus
vm_pu | va_degree | p_mw | q_mvar | lam_p | lam_q | |
---|---|---|---|---|---|---|
0 | 1.0 | 0.000000 | -23.521837 | NaN | 0.0 | 0.0 |
1 | 1.0 | -2.783612 | 60.000000 | 0.0 | 0.0 | 0.0 |
2 | 1.0 | 0.195199 | -29.310098 | NaN | 0.0 | 0.0 |
3 | 1.0 | -1.159257 | -7.168065 | NaN | 0.0 | 0.0 |
4 | 1.0 | -1.264625 | 0.000000 | 0.0 | 0.0 | 0.0 |
The julia file that is used to do that can be found in pandapower/pandapower/opf/run_powermodels_dc.jl and looks like this:
"""
using PowerModels
using Ipopt
using PP2PM
function run_powermodels(json_path)
pm = PP2PM.load_pm_from_json(json_path)
result = PowerModels.run_dc_opf(pm, Ipopt.IpoptSolver(),
setting = Dict("output" => Dict("branch_flows" => true)))
return result
end
""";
Of course PowerModels is a great modular tool that allows you to do much more than that. You might want to use a different OPF formulation, a relaxation method or a different solver. You might even want to use one of the variants of PowerModels that are being developed, such as PowerModelsACDC.jl or PowerModelsReliability.jl.
To do that, you can switch out the standard file with your own custom .jl file. Lets say we want to run a power flow instead of an OPF. There is a custom julia file for that in pandapower/tutorials/run_powermodels_custom.jl that looks like this:
"""
using PowerModels
using Ipopt
import JSON
function run_powermodels(json_path)
pm = PP2PM.load_pm_from_json(json_path)
result = PowerModels.run_pf(pm, ACPPowerModel, Ipopt.IpoptSolver(),
setting = Dict("output" => Dict("branch_flows" => true)))
return result
end
""";
We point the runpm function to this file, and as we can see by the flat voltage values, the OPF is now run with a DC network formulation:
pp.runpm(net, julia_file="run_powermodels_custom.jl")
net.res_bus
vm_pu | va_degree | p_mw | q_mvar | lam_p | lam_q | |
---|---|---|---|---|---|---|
0 | 1.010000 | -6.034381e-36 | -51.203353 | -1.245026 | 0.0 | 0.0 |
1 | 1.004756 | -5.426418e+00 | 60.000000 | 0.000000 | 0.0 | 0.0 |
2 | 1.010000 | -6.495320e+00 | 30.000000 | -15.224075 | 0.0 | 0.0 |
3 | 1.010000 | -3.621699e+00 | -40.000000 | 16.142383 | 0.0 | 0.0 |
4 | 1.010243 | -3.216964e+00 | 0.000000 | 0.000000 | 0.0 | 0.0 |
The PowerModels data structure that was passed to Julia can be accessed like this:
net._pm["bus"]
{'1': {'index': 1, 'bus_i': 1, 'zone': 1, 'bus_type': 3, 'vmax': 1.05, 'vmin': 0.95, 'va': 0.0, 'vm': 1.01, 'base_kv': 220.0}, '2': {'index': 2, 'bus_i': 2, 'zone': 1, 'bus_type': 1, 'vmax': 1.05, 'vmin': 0.95, 'va': 0.0, 'vm': 1.0, 'base_kv': 110.0}, '3': {'index': 3, 'bus_i': 3, 'zone': 1, 'bus_type': 2, 'vmax': 1.05, 'vmin': 0.95, 'va': 0.0, 'vm': 1.01, 'base_kv': 110.0}, '4': {'index': 4, 'bus_i': 4, 'zone': 1, 'bus_type': 2, 'vmax': 1.05, 'vmin': 0.95, 'va': 0.0, 'vm': 1.01, 'base_kv': 110.0}, '5': {'index': 5, 'bus_i': 5, 'zone': 1, 'bus_type': 1, 'vmax': 1.05, 'vmin': 0.95, 'va': 0.0, 'vm': 1.0, 'base_kv': 110.0}, '6': {'index': 6, 'bus_i': 6, 'zone': 1, 'bus_type': 1, 'vmax': 1.05, 'vmin': 0.95, 'va': 0.0, 'vm': 1.0, 'base_kv': 220.0}}
There is also a callback that allows you to add additional data to the PowerModels data structure in case it is not already added by the pandapower/PowerModels interface. In the callback you can add any data from the net, ppc or any source:
def add_data(net, ppc, pm):
pm["gen"]["1"]["bar"] = "foo"
pm["f_hz"] = net.f_hz
pp.runpm(net, julia_file="run_powermodels_custom.jl", pp_to_pm_callback=add_data)
print(net._pm["gen"]["1"]["bar"])
print(net._pm["f_hz"])
foo 50.0
These variables can now also be accessed on the Julia side, in case you need some more variables for custom optimizations.
Keep in mind that indices in PowerModels are 1-based so that the indices are shifted by one between ppc and pm. Furthermore, the net might contain some elements that are not in the ppc, as they are out of service or disconnected, which is why the indices of all elements have to be identified through the lookup tables in net._pd2ppc_lookups
.
Some notes on the internal data structure can be found in internal_datastructure.ipynb
Comparing the runopp function (that runs an OPF through PYPOWER) and the runpm function shows that PowerModels is much more performant:
%timeit pp.runopp(net)
1.22 s ± 14.9 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
%timeit pp.runpm_ac_opf(net)
221 ms ± 7.18 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)