Power flow#
SHOP implements PTDF-based power flow. The busbar and ac_line objects are used to specify the grid topology of the model. The following example demonstrates a simple model with power flow and thermal production units.
Set up model#
First, we initialize the model
import pandas as pd
import plotly.graph_objs as go
from pyshop import ShopSession
shop = ShopSession()
starttime=pd.Timestamp("2022-02-27T00")
endtime=starttime+pd.DateOffset(hours=24)
shop.set_time_resolution(starttime=starttime, endtime=endtime, timeunit="hour")
Add busbars and set the load. The load in busbar 1 is linearly increasing over time in this example, while all other input data is constant.
b1=shop.model.busbar.add_object("Busbar1")
b1.load.set(pd.Series([50+i for i in range(24)],index=[starttime+pd.Timedelta(hours=i) for i in range(24)]))
b2=shop.model.busbar.add_object("Busbar2")
b2.load.set(60)
b3=shop.model.busbar.add_object("Busbar3")
b3.load.set(300)
Then we define one thermal unit for each busbar and connect them. The marginal cost of the thermal units is equal to the fuel_cost in this case since n_segments is set to 1. Busbar 3 has the largest thermal unit, while busbar 1 has the cheapest one. Only a small and expensive unit is located in busbar 2.
t1=shop.model.thermal.add_object("Thermal1")
t1.fuel_cost.set(75)
t1.n_segments.set(1)
t1.min_prod.set(0)
t1.max_prod.set(140)
t1.connect_to(b1)
t2=shop.model.thermal.add_object("Thermal2")
t2.fuel_cost.set(180)
t2.n_segments.set(1)
t2.min_prod.set(0)
t2.max_prod.set(20)
t2.connect_to(b2)
t3=shop.model.thermal.add_object("Thermal3")
t3.fuel_cost.set(140)
t3.n_segments.set(1)
t3.min_prod.set(0)
t3.max_prod.set(400)
t3.connect_to(b3)
There are three AC lines that connect the busbars, and there is no power loss in the transmission. The AC line connecting busbar 2 and 3 has the lowest transmission capacity, which is specified in both directions by the max_forward_flow and max_backward_flow attributes. Note that the positive (forward) flow direction on each line is determined by the connections between the busbar and ac_line objects.
Create AC lines and connect between buses:
ac1=shop.model.ac_line.add_object("AC_line1")
ac1.max_forward_flow.set(40)
ac1.max_backward_flow.set(40)
b1.connect_to(ac1)
ac1.connect_to(b2)
ac2=shop.model.ac_line.add_object("AC_line2")
ac2.max_forward_flow.set(250)
ac2.max_backward_flow.set(250)
b1.connect_to(ac2)
ac2.connect_to(b3)
ac3=shop.model.ac_line.add_object("AC_line3")
ac3.max_forward_flow.set(10)
ac3.max_backward_flow.set(10)
b2.connect_to(ac3)
ac3.connect_to(b3)
The we set the PTDF matrix that determines the flow on each line as a function of the net power injection in each busbar. Busbar 3 is chosen as the slack bus in this example, which means that the PTDF values for busbar 3 are all zero.
ac_lines=["AC_line1", "AC_line2", "AC_line3"]
b1.ptdf.set(pd.Series(data=[0.2/0.5, (0.2+0.1)/0.5, 0.2/0.5], index=ac_lines))
b2.ptdf.set(pd.Series(data=[-0.1/0.5, 0.1/0.5, (0.2+0.2)/0.5], index=ac_lines))
b3.ptdf.set(pd.Series(data=[0.0, 0.0, 0.0], index=ac_lines))
The grid topology with the PTDFs visualized for busbars 1 and 2 are shown in the figures below.
|
|---|
Grid topology with the PTDF for busbar 1 shown as red numbers. |
|
|---|
Grid topology with the PTDF for busbar 2 shown as green numbers. Note that the PTDF for line 2 is negative due to the chosen positive line direction. |
Run the model:
shop.start_sim([],["1"])
Results#
First, we take a look at the production of each unit in the system. The cheapest unit located in busbar 2 does not operate at full capacity until 15 hours after optimization start, this is due to the transmission line constraints and PTDF values used in the example.
p1 = t1.production.get()
p2 = t2.production.get()
p3 = t3.production.get()
fig = go.Figure()
fig.add_trace(go.Scatter(x=p1.index, y=p1.values, name="Sum production busbar 1"))
fig.add_trace(go.Scatter(x=p2.index, y=p2.values, name="Sum production busbar 2"))
fig.add_trace(go.Scatter(x=p3.index, y=p3.values, name="Sum production busbar 3"))
fig.update_layout(xaxis_title="Time", yaxis_title="Production [MW]")
fig.show()
/opt/conda/lib/python3.13/site-packages/kaleido/_sync_server.py:11: UserWarning:
Warning: You have Plotly version 5.24.0, which is not compatible with this version of Kaleido (1.5.0).
This means that static image generation (e.g. `fig.write_image()`) will not work.
Please upgrade Plotly to version 6.1.1 or greater, or downgrade Kaleido to version 0.2.1.
You can however, use the Kaleido API directly which will work with your plotly version. `kaleido.write_fig(...)`, for example. Please see the kaleido documentation.
The total net power injection in each busbar determines the flow on each line in combination with the PTDF matrix.
l1 = b1.load.get()
l2 = b2.load.get()
l3 = b3.load.get()
fig = go.Figure()
fig.add_trace(go.Scatter(x=l1.index, y=p1.values-l1.values, name="Net injection busbar 1"))
fig.add_trace(go.Scatter(x=l2.index, y=p2.values-l2.values, name="Net injection busbar 2"))
fig.add_trace(go.Scatter(x=l3.index, y=p3.values-l3.values, name="Net injection busbar 3"))
fig.update_layout(xaxis_title="Time", yaxis_title="Power [MW]")
fig.show()
The resulting flow on each power line is shown below. Note that the line between busbar 2 and 3 operates at full (reverse) capacity for all times.
flow_12 = ac1.flow.get()
flow_13 = ac2.flow.get()
flow_23 = ac3.flow.get()
fig = go.Figure()
fig.add_trace(go.Scatter(x=flow_12.index, y=flow_12.values, name="Flow from 1 to 2"))
fig.add_trace(go.Scatter(x=flow_13.index, y=flow_13.values, name="Flow from 1 to 3"))
fig.add_trace(go.Scatter(x=flow_23.index, y=flow_23.values, name="Flow from 2 to 3"))
fig.update_layout(xaxis_title="Time", yaxis_title="Power flow [MW]")
fig.show()
Finally, we can plot the dual value of the power balance in each busbar by retrieving the energy_price attribute, which represents the nodal price. Busbars 2 and 3 have a constant price equal to the marginal cost of their respective thermal units. This is also true for busbar 1 initially when its thermal unit operates below maximum capacity. From hour 15, the price in busbar 1 jumps up to 160 €/MWh, which is the average of the price in the other busbars. This is caused by the constrained line between busbar 2 and 3, which forces the increasing load in busbar 1 to be served by increasing the production of both the thermal unit in busbar 2 and 3 by the same amount. The marginal cost of increasing the load in busbar 1 is therefore 50% of the marginal cost of unit 2 plus 50% of the marginal cost of unit 3.
price_1 = b1.energy_price.get()
price_2 = b2.energy_price.get()
price_3 = b3.energy_price.get()
fig = go.Figure()
fig.add_trace(go.Scatter(x=price_1.index, y=price_1.values, name="Price in busbar 1"))
fig.add_trace(go.Scatter(x=price_2.index, y=price_2.values, name="Price in busbar 2"))
fig.add_trace(go.Scatter(x=price_3.index, y=price_3.values, name="Price in busbar 3"))
fig.update_layout(xaxis_title="Time", yaxis_title="Power price [€/MWh]")
fig.show()