Security-Constrained LOPF¶
In this example, the dispatch of generators is optimised using the security-constrained linear optimal power flow (SCLOPF) functionality of PyPSA, to guarantee that no branches are overloaded in the event of certain branch outages.
In [2]:
Copied!
import pypsa
n = pypsa.examples.scigrid_de()
# correct some infeasibilties in the network
for line_name in ["316", "527", "602"]:
n.lines.loc[line_name, "s_nom"] = 1200
now = n.snapshots[0]
n.plot(bus_size=0);
import pypsa
n = pypsa.examples.scigrid_de()
# correct some infeasibilties in the network
for line_name in ["316", "527", "602"]:
n.lines.loc[line_name, "s_nom"] = 1200
now = n.snapshots[0]
n.plot(bus_size=0);
INFO:pypsa.network.io:Imported network 'SciGrid-DE' has buses, carriers, generators, lines, loads, storage_units, transformers
First, let's run the network without any $N-1$ security constraints and see how much it costs to operate the system.:
In [3]:
Copied!
n0 = n.copy()
n0.optimize(snapshots=now)
n0.statistics.opex().sum()
n0 = n.copy()
n0.optimize(snapshots=now)
n0.statistics.opex().sum()
/tmp/ipykernel_8061/344066763.py:2: 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.
n0.optimize(snapshots=now)
WARNING:pypsa.consistency:The following transformers have zero r, which could break the linear load flow:
Index(['2', '5', '10', '12', '13', '15', '18', '20', '22', '24', '26', '30',
'32', '37', '42', '46', '52', '56', '61', '68', '69', '74', '78', '86',
'87', '94', '95', '96', '99', '100', '104', '105', '106', '107', '117',
'120', '123', '124', '125', '128', '129', '138', '143', '156', '157',
'159', '160', '165', '184', '191', '195', '201', '220', '231', '232',
'233', '236', '247', '248', '250', '251', '252', '261', '263', '264',
'267', '272', '279', '281', '282', '292', '303', '307', '308', '312',
'315', '317', '322', '332', '334', '336', '338', '351', '353', '360',
'362', '382', '384', '385', '391', '403', '404', '413', '421', '450',
'458'],
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.08s
INFO:linopy.constants: Optimization successful: Status: ok Termination condition: optimal Solution: 2485 primals, 5957 duals Objective: 3.32e+05 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, Transformer-fix-s-lower, Transformer-fix-s-upper, StorageUnit-fix-p_dispatch-lower, StorageUnit-fix-p_dispatch-upper, StorageUnit-fix-p_store-lower, StorageUnit-fix-p_store-upper, StorageUnit-fix-state_of_charge-lower, StorageUnit-fix-state_of_charge-upper, Kirchhoff-Voltage-Law, StorageUnit-energy_balance were not assigned to the network.
Out[3]:
np.float64(332383.51399999997)
The security-constrained linear optimal power flow (SCLOPF) is executed using the n.optimize.optimize_security_constrained() method.
This method takes a list of snapshots (here: now) and a list of branches that are to be considered as outages (here, the 30 lines with the highest loading).
In [4]:
Copied!
branch_outages = (n0.lines_t.p0.loc[now] / n0.lines.s_nom).nlargest(30).index
n.optimize.optimize_security_constrained(
now, branch_outages=branch_outages, log_to_console=False
)
n.statistics.opex().sum()
branch_outages = (n0.lines_t.p0.loc[now] / n0.lines.s_nom).nlargest(30).index
n.optimize.optimize_security_constrained(
now, branch_outages=branch_outages, log_to_console=False
)
n.statistics.opex().sum()
WARNING:pypsa.consistency:The following transformers have zero r, which could break the linear load flow:
Index(['2', '5', '10', '12', '13', '15', '18', '20', '22', '24', '26', '30',
'32', '37', '42', '46', '52', '56', '61', '68', '69', '74', '78', '86',
'87', '94', '95', '96', '99', '100', '104', '105', '106', '107', '117',
'120', '123', '124', '125', '128', '129', '138', '143', '156', '157',
'159', '160', '165', '184', '191', '195', '201', '220', '231', '232',
'233', '236', '247', '248', '250', '251', '252', '261', '263', '264',
'267', '272', '279', '281', '282', '292', '303', '307', '308', '312',
'315', '317', '322', '332', '334', '336', '338', '351', '353', '360',
'362', '382', '384', '385', '391', '403', '404', '413', '421', '450',
'458'],
dtype='object', name='name')
/home/docs/checkouts/readthedocs.org/user_builds/pypsa/checkouts/latest/pypsa/optimization/abstract.py:436: 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. m = n.optimize.create_model(
INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.model:Solver options: - log_to_console: False
INFO:linopy.io: Writing time: 0.13s
INFO:linopy.constants: Optimization successful: Status: ok Termination condition: optimal Solution: 2485 primals, 62837 duals Objective: 4.25e+05 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, Transformer-fix-s-lower, Transformer-fix-s-upper, StorageUnit-fix-p_dispatch-lower, StorageUnit-fix-p_dispatch-upper, StorageUnit-fix-p_store-lower, StorageUnit-fix-p_store-upper, StorageUnit-fix-state_of_charge-lower, StorageUnit-fix-state_of_charge-upper, Kirchhoff-Voltage-Law, StorageUnit-energy_balance, Line-fix-s-lower-security-for-Line-outage-in-sub-network-0, Line-fix-s-upper-security-for-Line-outage-in-sub-network-0, Transformer-fix-s-lower-security-for-Line-outage-in-sub-network-0, Transformer-fix-s-upper-security-for-Line-outage-in-sub-network-0 were not assigned to the network.
Out[4]:
np.float64(424597.83679)
You can see that the cost of operating the system in the given hour rises from 332 k€ to 427 k€ when the security constraints are applied.
The maps below indicate the difference in the dispatch patterns, line flows and marginal prices between the $N-0$ (first) and $N-1$ (second) cases.
In [5]:
Copied!
def plot_network(n, snapshot):
bus_size = (
n.statistics.supply(groupby="bus", components=["Generator", "StorageUnit"])
.groupby("bus")
.sum()
)
line_flows = n.lines_t.p0.loc[snapshot]
bus_color = n.buses_t.marginal_price.loc[snapshot]
line_loading = n.lines_t.p0.abs().loc[snapshot] / n.lines.s_nom
n.plot(
bus_size=bus_size / 30000,
bus_color=bus_color,
bus_cmap="Reds",
line_color=line_loading,
line_flow=line_flows / 50,
)
def plot_network(n, snapshot):
bus_size = (
n.statistics.supply(groupby="bus", components=["Generator", "StorageUnit"])
.groupby("bus")
.sum()
)
line_flows = n.lines_t.p0.loc[snapshot]
bus_color = n.buses_t.marginal_price.loc[snapshot]
line_loading = n.lines_t.p0.abs().loc[snapshot] / n.lines.s_nom
n.plot(
bus_size=bus_size / 30000,
bus_color=bus_color,
bus_cmap="Reds",
line_color=line_loading,
line_flow=line_flows / 50,
)
In [6]:
Copied!
plot_network(n0, now)
plot_network(n0, now)
In [7]:
Copied!
plot_network(n, now)
plot_network(n, now)
We can also look at where the nodal dispatch is ramped up (red) or down (blue) in the $N-1$ case compared to the $N-0$ case. Mostly this means ramp down upstream of the potential outages, and ramp up downstream of the potential outages.
In [8]:
Copied!
bus_size0 = (
(
n0.statistics.supply(groupby="bus", components=["Generator", "StorageUnit"])
.groupby("bus")
.sum()
)
.reindex(index=n.buses.index)
.fillna(0)
)
bus_size1 = (
(
n.statistics.supply(groupby="bus", components=["Generator", "StorageUnit"])
.groupby("bus")
.sum()
)
.reindex(index=n.buses.index)
.fillna(0)
)
bus_size = bus_size1 - bus_size0
n.plot(
bus_size=bus_size.abs() / 30000,
bus_color=bus_size.map(lambda x: "b" if x < 0 else "r"),
);
bus_size0 = (
(
n0.statistics.supply(groupby="bus", components=["Generator", "StorageUnit"])
.groupby("bus")
.sum()
)
.reindex(index=n.buses.index)
.fillna(0)
)
bus_size1 = (
(
n.statistics.supply(groupby="bus", components=["Generator", "StorageUnit"])
.groupby("bus")
.sum()
)
.reindex(index=n.buses.index)
.fillna(0)
)
bus_size = bus_size1 - bus_size0
n.plot(
bus_size=bus_size.abs() / 30000,
bus_color=bus_size.map(lambda x: "b" if x < 0 else "r"),
);
We can also double-check that the $N-1$ constraints are satisfied and no lines are overloaded in the event of the outages considered:
In [9]:
Copied!
n.optimize.fix_optimal_dispatch()
p0_test = n.lpf_contingency(now, branch_outages=branch_outages)
n.optimize.fix_optimal_dispatch()
p0_test = n.lpf_contingency(now, branch_outages=branch_outages)
INFO:pypsa.network.power_flow:Performing linear load-flow on AC sub-network <pypsa.SubNetwork object at 0x7f8c951396d0> for snapshot(s) DatetimeIndex(['2011-01-01'], dtype='datetime64[ns]', name='snapshot', freq=None)
WARNING:pypsa.network.power_flow:No type given for 183, assuming it is a line
WARNING:pypsa.network.power_flow:No type given for 214, assuming it is a line
WARNING:pypsa.network.power_flow:No type given for 350, assuming it is a line
WARNING:pypsa.network.power_flow:No type given for 389, assuming it is a line
WARNING:pypsa.network.power_flow:No type given for 448, assuming it is a line
WARNING:pypsa.network.power_flow:No type given for 586, assuming it is a line
WARNING:pypsa.network.power_flow:No type given for 608, assuming it is a line
WARNING:pypsa.network.power_flow:No type given for 707, assuming it is a line
WARNING:pypsa.network.power_flow:No type given for 809, assuming it is a line
WARNING:pypsa.network.power_flow:No type given for 818, assuming it is a line
WARNING:pypsa.network.power_flow:No type given for 390, assuming it is a line
WARNING:pypsa.network.power_flow:No type given for 587, assuming it is a line
WARNING:pypsa.network.power_flow:No type given for 124, assuming it is a line
WARNING:pypsa.network.power_flow:No type given for 682, assuming it is a line
WARNING:pypsa.network.power_flow:No type given for 683, assuming it is a line
WARNING:pypsa.network.power_flow:No type given for 29, assuming it is a line
WARNING:pypsa.network.power_flow:No type given for 218, assuming it is a line
WARNING:pypsa.network.power_flow:No type given for 714, assuming it is a line
WARNING:pypsa.network.power_flow:No type given for 595, assuming it is a line
WARNING:pypsa.network.power_flow:No type given for 596, assuming it is a line
WARNING:pypsa.network.power_flow:No type given for 371, assuming it is a line
WARNING:pypsa.network.power_flow:No type given for 17, assuming it is a line
WARNING:pypsa.network.power_flow:No type given for 514, assuming it is a line
WARNING:pypsa.network.power_flow:No type given for 511, assuming it is a line
WARNING:pypsa.network.power_flow:No type given for 510, assuming it is a line
WARNING:pypsa.network.power_flow:No type given for 58, assuming it is a line
WARNING:pypsa.network.power_flow:No type given for 254, assuming it is a line
WARNING:pypsa.network.power_flow:No type given for 792, assuming it is a line
WARNING:pypsa.network.power_flow:No type given for 436, assuming it is a line
WARNING:pypsa.network.power_flow:No type given for 11, assuming it is a line
Check the maximum loading as per unit of s_nom in each contingency:
In [10]:
Copied!
max_loading = abs(p0_test.divide(n.passive_branches().s_nom, axis=0)).max()
max_loading
max_loading = abs(p0_test.divide(n.passive_branches().s_nom, axis=0)).max()
max_loading
Out[10]:
base 1.0 (Line, 183) 1.0 (Line, 214) 1.0 (Line, 350) 1.0 (Line, 389) 1.0 (Line, 448) 1.0 (Line, 586) 1.0 (Line, 608) 1.0 (Line, 707) 1.0 (Line, 809) 1.0 (Line, 818) 1.0 (Line, 390) 1.0 (Line, 587) 1.0 (Line, 124) 0.0 (Line, 682) 1.0 (Line, 683) 1.0 (Line, 29) 1.0 (Line, 218) 1.0 (Line, 714) 1.0 (Line, 595) 1.0 (Line, 596) 1.0 (Line, 371) 1.0 (Line, 17) 1.0 (Line, 514) 1.0 (Line, 511) 1.0 (Line, 510) 1.0 (Line, 58) 1.0 (Line, 254) 1.0 (Line, 792) 1.0 (Line, 436) 1.0 (Line, 11) 1.0 dtype: float64
In [ ]:
Copied!