Quickstart 2 - Power Flow¶
Problem Description¶
Consider a three-zone electricity market with the following loads and generators:
| Zone | Load (MW) | Generator | Capacity (MW) | Marginal Cost (€/MWh) |
|---|---|---|---|---|
| 1 | 50 | Gen A | 140 | 7.5 |
| Gen B | 285 | 6.0 | ||
| 2 | 60 | Gen C | 90 | 14.0 |
| 3 | 300 | Gen D | 85 | 10.0 |
The zones are connected by transmission lines running on 11 kV with the following characteristics:
| Line | From | To | Capacity (MW) | Reactance (Ohm) |
|---|---|---|---|---|
| 1 | Zone 1 | Zone 2 | 126 | 0.2 |
| 2 | Zone 1 | Zone 3 | 250 | 0.2 |
| 3 | Zone 2 | Zone 3 | 130 | 0.1 |
Find the least-cost dispatch of the generators to meet the loads while respecting the transmission capacities. Identify the marginal prices at each bus and the flow on the transmission lines. Calculate the non-linear power flow based on the optimised results, and determine the losses on each line as well as the voltage angles at each bus.
PyPSA Solution¶
Start by importing PyPSA, creating a new network instance and adding the three buses with voltage level 11 kV:
from numpy import pi
import pypsa
n = pypsa.Network()
n.add("Bus", ["zone_1", "zone_2", "zone_3"], v_nom=11)
Adding multiple components with n.add() is supported by passing lists for each argument. Any scalar values are broadcasted to all components.
Thus, we can add the buses, loads, generators, and lines in one go and solve the network::
n.add(
"Load",
["load_1", "load_2", "load_3"],
bus=["zone_1", "zone_2", "zone_3"],
p_set=[50, 60, 300],
)
n.add(
"Generator",
["gen_A", "gen_B", "gen_C", "gen_D"],
bus=["zone_1", "zone_1", "zone_2", "zone_3"],
p_nom=[140, 285, 90, 85],
marginal_cost=[7.5, 6, 14, 10],
)
n.add(
"Line",
["line_1", "line_2", "line_3"],
bus0=["zone_1", "zone_1", "zone_2"],
bus1=["zone_2", "zone_3", "zone_3"],
s_nom=[126, 250, 130],
x=[0.02, 0.02, 0.01],
r=0.01,
)
n.optimize()
/tmp/ipykernel_3829/820350409.py:26: FutureWarning: The default value of `include_objective_constant` will change from True to False in version 2.0. Set `include_objective_constant` explicitly to suppress this warning. Using False improves LP numerical conditioning by not including the objective constant as a variable. n.optimize() WARNING:pypsa.consistency:The following buses have carriers which are not defined. Run n.sanitize() to add them. Components with undefined carriers: Index(['zone_1', 'zone_2', 'zone_3'], dtype='object', name='name')
WARNING:pypsa.consistency:The following lines have carriers which are not defined. Run n.sanitize() to add them. Components with undefined carriers: Index(['line_1', 'line_2', 'line_3'], dtype='object', name='name')
INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.model:Solver options: - log_to_console: False
INFO:linopy.io: Writing time: 0.03s
INFO:linopy.constants: Optimization successful: Status: ok Termination condition: optimal Solution: 7 primals, 18 duals Objective: 2.84e+03 Solver model: available Solver message: Optimal
INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-fix-p-lower, Generator-fix-p-upper, Line-fix-s-lower, Line-fix-s-upper, Kirchhoff-Voltage-Law were not assigned to the network.
('ok', 'optimal')
It is possible to inspect the underlying linopy model once either n.optimize() or n.optimize.create_model() have been called. This can be done by accessing the n.model attribute, which prints the variable and constraint names and their dimensions.
n.model
Linopy LP model =============== Variables: ---------- * Generator-p (snapshot, name) * Line-s (snapshot, name) Constraints: ------------ * Generator-fix-p-lower (snapshot, name) * Generator-fix-p-upper (snapshot, name) * Line-fix-s-lower (snapshot, name) * Line-fix-s-upper (snapshot, name) * Bus-nodal_balance (name, snapshot) * Kirchhoff-Voltage-Law (snapshot, cycle) Status: ------- ok
To look at the equations of a specific constraint, such as the nodal balances of the buses in each time step, execute:
n.model.constraints["Bus-nodal_balance"]
Constraint `Bus-nodal_balance` [name: 3, snapshot: 1]: ------------------------------------------------------ [zone_1, now]: +1 Generator-p[now, gen_A] + 1 Generator-p[now, gen_B] - 1 Line-s[now, line_1] - 1 Line-s[now, line_2] = 50.0 [zone_2, now]: +1 Generator-p[now, gen_C] - 1 Line-s[now, line_3] + 1 Line-s[now, line_1] = 60.0 [zone_3, now]: +1 Generator-p[now, gen_D] + 1 Line-s[now, line_2] + 1 Line-s[now, line_3] = 300.0
As in the first example, the optimal solution can be accessed from the network object:
display(n.buses_t.marginal_price)
display(n.generators_t.p)
display(n.lines_t.p0)
| name | zone_1 | zone_2 | zone_3 |
|---|---|---|---|
| snapshot | |||
| now | 7.5 | 11.25 | 10.0 |
| name | gen_A | gen_B | gen_C | gen_D |
|---|---|---|---|---|
| snapshot | ||||
| now | 50.0 | 285.0 | -0.0 | 75.0 |
| name | line_1 | line_2 | line_3 |
|---|---|---|---|
| snapshot | |||
| now | 126.0 | 159.0 | 66.0 |
For the next step, we want to calculate the non-linear AC power flow based on the optimised results (which uses the linearised DC power flow model) using the Newton-Raphson method.
This can be done by calling n.pf(), but first we need to provide the set points for the generator dispatch.
These will be kept in the power flow, except for the slack generator, the output of which will be adjusted to balance the network.
n.optimize.fix_optimal_dispatch()
display(n.generators_t.p_set)
display(n.generators.control)
n.pf()
| name | gen_A | gen_B | gen_C | gen_D |
|---|---|---|---|---|
| snapshot | ||||
| now | 50.0 | 285.0 | -0.0 | 75.0 |
name gen_A Slack gen_B PQ gen_C PQ gen_D PQ Name: control, dtype: object
INFO:pypsa.network.power_flow:Performing non-linear load-flow on AC sub-network <pypsa.SubNetwork object at 0x7d77c71f6ad0> for snapshots Index(['now'], dtype='object', name='snapshot')
{'n_iter': name 0
snapshot
now 3,
'error': name 0
snapshot
now 1.491031e-09,
'converged': name 0
snapshot
now True}
Alright, it converged! Now we can see how the slack generator is used to balance the line losses as well as the reactive power of the generators and lines. The linear power flow necglects both.
display(n.generators_t.p)
display(n.generators_t.q)
display(n.lines_t.q0)
display((n.lines_t.p0 + n.lines_t.p1).sum().sum()) # active power losses
| name | gen_A | gen_B | gen_C | gen_D |
|---|---|---|---|---|
| snapshot | ||||
| now | 53.860705 | 285.0 | -0.0 | 75.0 |
| name | gen_A | gen_B | gen_C | gen_D |
|---|---|---|---|---|
| snapshot | ||||
| now | 7.379326 | 0.0 | 0.0 | 0.0 |
| name | line_1 | line_2 | line_3 |
|---|---|---|---|
| snapshot | |||
| now | -1.811567 | 9.190894 | -4.388279 |
np.float64(3.86070466597905)
You can see that the increase in the slack generator's output is equal to the sum of the active power losses in the lines.
With non-linear power flow, the voltage angles at the buses are calculated as well, which are stored in radians but can easily be converted to degrees:
n.buses_t.v_ang * 180 / pi
| name | zone_1 | zone_2 | zone_3 |
|---|---|---|---|
| snapshot | |||
| now | 0.0 | -1.202765 | -1.532529 |
This example is based on Tom Brown's Energy Systems course, taken from the lecture Complex Markets, slides 42ff.