Note
You can download this example as a Jupyter notebook or start it in interactive mode.
Optimization without pyomo
In this example we demonstrate the behaviour of the Linear Optimal Power Flow (LOPF) calculation without using pyomo. This requires to set pyomo to False in the lopf function. Then, the communication with the solvers happens via in-house functions which leads to a much faster solving process.
[1]:
import pypsa
import pandas as pd
import os
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
Cell In [1], line 1
----> 1 import pypsa
2 import pandas as pd
3 import os
File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.21.2/lib/python3.11/site-packages/pypsa/__init__.py:10
1 # -*- coding: utf-8 -*-
4 """
5 Python for Power Systems Analysis (PyPSA)
6
7 Grid calculation library.
8 """
---> 10 from pypsa import (
11 components,
12 contingency,
13 descriptors,
14 examples,
15 geo,
16 io,
17 linopf,
18 linopt,
19 networkclustering,
20 opf,
21 opt,
22 optimization,
23 pf,
24 plot,
25 )
26 from pypsa.components import Network, SubNetwork
28 __version__ = "0.21.2"
File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.21.2/lib/python3.11/site-packages/pypsa/components.py:50
37 from pypsa.io import (
38 export_to_csv_folder,
39 export_to_hdf5,
(...)
47 import_series_from_dataframe,
48 )
49 from pypsa.opf import network_lopf, network_opf
---> 50 from pypsa.optimization.optimize import OptimizationAccessor
51 from pypsa.pf import (
52 calculate_B_H,
53 calculate_dependent_values,
(...)
62 sub_network_pf,
63 )
64 from pypsa.plot import iplot, plot
File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.21.2/lib/python3.11/site-packages/pypsa/optimization/__init__.py:7
1 #!/usr/bin/env python3
2 # -*- coding: utf-8 -*-
3 """
4 Build optimisation problems from PyPSA networks with Linopy.
5 """
----> 7 from pypsa.optimization import abstract, constraints, optimize, variables
8 from pypsa.optimization.optimize import create_model
File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.21.2/lib/python3.11/site-packages/pypsa/optimization/constraints.py:9
6 import logging
8 import pandas as pd
----> 9 from linopy.expressions import LinearExpression, merge
10 from numpy import arange, cumsum, inf, nan, roll
11 from scipy import sparse
File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.21.2/lib/python3.11/site-packages/linopy/__init__.py:9
1 #!/usr/bin/env python3
2 # -*- coding: utf-8 -*-
3 """
4 Created on Wed Mar 10 11:03:06 2021.
5
6 @author: fabulous
7 """
----> 9 from linopy import model, remote
10 from linopy.expressions import merge
11 from linopy.io import read_netcdf
File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.21.2/lib/python3.11/site-packages/linopy/model.py:22
20 from linopy import solvers
21 from linopy.common import best_int, replace_by_map
---> 22 from linopy.constraints import (
23 AnonymousConstraint,
24 AnonymousScalarConstraint,
25 Constraints,
26 )
27 from linopy.eval import Expr
28 from linopy.expressions import LinearExpression, ScalarLinearExpression
File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.21.2/lib/python3.11/site-packages/linopy/constraints.py:21
18 from scipy.sparse import coo_matrix
19 from xarray import DataArray, Dataset
---> 21 from linopy import expressions, variables
22 from linopy.common import (
23 _merge_inplace,
24 has_assigned_model,
(...)
27 replace_by_map,
28 )
31 class Constraint(DataArray):
File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.21.2/lib/python3.11/site-packages/linopy/expressions.py:23
20 from xarray.core.dataarray import DataArrayCoordinates
21 from xarray.core.groupby import _maybe_reorder, peek_at
---> 23 from linopy import constraints, variables
24 from linopy.common import as_dataarray
27 def exprwrap(method, *default_args, **new_default_kwargs):
File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.21.2/lib/python3.11/site-packages/linopy/variables.py:398
393 roll = varwrap(DataArray.roll)
395 rolling = varwrap(DataArray.rolling)
--> 398 @dataclass(repr=False)
399 class Variables:
400 """
401 A variables container used for storing multiple variable arrays.
402 """
404 labels: Dataset = Dataset()
File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.21.2/lib/python3.11/dataclasses.py:1211, in dataclass.<locals>.wrap(cls)
1210 def wrap(cls):
-> 1211 return _process_class(cls, init, repr, eq, order, unsafe_hash,
1212 frozen, match_args, kw_only, slots,
1213 weakref_slot)
File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.21.2/lib/python3.11/dataclasses.py:959, in _process_class(cls, init, repr, eq, order, unsafe_hash, frozen, match_args, kw_only, slots, weakref_slot)
956 kw_only = True
957 else:
958 # Otherwise it's a field of some type.
--> 959 cls_fields.append(_get_field(cls, name, type, kw_only))
961 for f in cls_fields:
962 fields[f.name] = f
File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.21.2/lib/python3.11/dataclasses.py:816, in _get_field(cls, a_name, a_type, default_kw_only)
812 # For real fields, disallow mutable defaults. Use unhashable as a proxy
813 # indicator for mutability. Read the __hash__ attribute from the class,
814 # not the instance.
815 if f._field_type is _FIELD and f.default.__class__.__hash__ is None:
--> 816 raise ValueError(f'mutable default {type(f.default)} for field '
817 f'{f.name} is not allowed: use default_factory')
819 return f
ValueError: mutable default <class 'xarray.core.dataset.Dataset'> for field labels is not allowed: use default_factory
[2]:
n = pypsa.examples.ac_dc_meshed(from_master=True)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In [2], line 1
----> 1 n = pypsa.examples.ac_dc_meshed(from_master=True)
NameError: name 'pypsa' is not defined
Modify the network a bit: We set gas generators to non-extendable
[3]:
n.generators.loc[n.generators.carrier == "gas", "p_nom_extendable"] = False
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In [3], line 1
----> 1 n.generators.loc[n.generators.carrier == "gas", "p_nom_extendable"] = False
NameError: name 'n' is not defined
Add ramp limit
[4]:
n.generators.loc[n.generators.carrier == "gas", "ramp_limit_down"] = 0.2
n.generators.loc[n.generators.carrier == "gas", "ramp_limit_up"] = 0.2
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In [4], line 1
----> 1 n.generators.loc[n.generators.carrier == "gas", "ramp_limit_down"] = 0.2
2 n.generators.loc[n.generators.carrier == "gas", "ramp_limit_up"] = 0.2
NameError: name 'n' is not defined
Add additional storage units (cyclic and non-cyclic) and fix one state_of_charge
[5]:
n.add(
"StorageUnit",
"su",
bus="Manchester",
marginal_cost=10,
inflow=50,
p_nom_extendable=True,
capital_cost=10,
p_nom=2000,
efficiency_dispatch=0.5,
cyclic_state_of_charge=True,
state_of_charge_initial=1000,
)
n.add(
"StorageUnit",
"su2",
bus="Manchester",
marginal_cost=10,
p_nom_extendable=True,
capital_cost=50,
p_nom=2000,
efficiency_dispatch=0.5,
carrier="gas",
cyclic_state_of_charge=False,
state_of_charge_initial=1000,
)
n.storage_units_t.state_of_charge_set.loc[n.snapshots[7], "su"] = 100
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In [5], line 1
----> 1 n.add(
2 "StorageUnit",
3 "su",
4 bus="Manchester",
5 marginal_cost=10,
6 inflow=50,
7 p_nom_extendable=True,
8 capital_cost=10,
9 p_nom=2000,
10 efficiency_dispatch=0.5,
11 cyclic_state_of_charge=True,
12 state_of_charge_initial=1000,
13 )
15 n.add(
16 "StorageUnit",
17 "su2",
(...)
26 state_of_charge_initial=1000,
27 )
29 n.storage_units_t.state_of_charge_set.loc[n.snapshots[7], "su"] = 100
NameError: name 'n' is not defined
Add an additional store
[6]:
n.add("Bus", "storebus", carrier="hydro", x=-5, y=55)
n.madd(
"Link",
["battery_power", "battery_discharge"],
"",
bus0=["Manchester", "storebus"],
bus1=["storebus", "Manchester"],
p_nom=100,
efficiency=0.9,
p_nom_extendable=True,
p_nom_max=1000,
)
n.madd(
"Store",
["store"],
bus="storebus",
e_nom=2000,
e_nom_extendable=True,
marginal_cost=10,
capital_cost=10,
e_nom_max=5000,
e_initial=100,
e_cyclic=True,
);
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In [6], line 1
----> 1 n.add("Bus", "storebus", carrier="hydro", x=-5, y=55)
2 n.madd(
3 "Link",
4 ["battery_power", "battery_discharge"],
(...)
11 p_nom_max=1000,
12 )
13 n.madd(
14 "Store",
15 ["store"],
(...)
23 e_cyclic=True,
24 );
NameError: name 'n' is not defined
Extra functionalities:
[7]:
from pypsa.linopt import get_var, linexpr, join_exprs, define_constraints
One of the most important functions is linexpr which take one or more tuples of coefficient and variable pairs which should go into the left hand side (lhs) of the constraint.
Add mimimum for state_of_charge
[8]:
def minimal_state_of_charge(n, snapshots):
vars_soc = get_var(n, "StorageUnit", "state_of_charge")
lhs = linexpr((1, vars_soc))
define_constraints(n, lhs, ">", 50, "StorageUnit", "soc_lower_bound")
Fix the ratio between ingoing and outgoing capacity of the Store
[9]:
def fix_link_cap_ratio(n, snapshots):
vars_link = get_var(n, "Link", "p_nom")
eff = n.links.at["battery_power", "efficiency"]
lhs = linexpr(
(1, vars_link["battery_power"]), (-eff, vars_link["battery_discharge"])
)
define_constraints(n, lhs, "=", 0, "battery_discharge", attr="fixratio")
Every bus must in total produce the 20% of the total demand
This requires the function pypsa.linopt.join_exprs which sums up arrays of linear expressions
[10]:
def fix_bus_production(n, snapshots):
total_demand = n.loads_t.p_set.sum().sum()
prod_per_bus = (
linexpr((1, get_var(n, "Generator", "p")))
.groupby(n.generators.bus, axis=1)
.apply(join_exprs)
)
define_constraints(
n, prod_per_bus, ">=", total_demand / 5, "Bus", "production_share"
)
Combine them …
[11]:
def extra_functionalities(n, snapshots):
minimal_state_of_charge(n, snapshots)
fix_link_cap_ratio(n, snapshots)
fix_bus_production(n, snapshots)
…and run the lopf with pyomo=False
[12]:
n.lopf(
pyomo=False,
extra_functionality=extra_functionalities,
keep_shadowprices=["Bus", "battery_discharge", "StorageUnit"],
)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In [12], line 1
----> 1 n.lopf(
2 pyomo=False,
3 extra_functionality=extra_functionalities,
4 keep_shadowprices=["Bus", "battery_discharge", "StorageUnit"],
5 )
NameError: name 'n' is not defined
The keep_shadowprices argument in the lopf now decides which shadow prices (SP) should be retrieved. It can either be set to True, then all SP are kept. It also can be a list of names of the constraints. Therefore the name argument in define_constraints is necessary, in our case ‘battery_discharge’, ‘StorageUnit’ and ‘Bus’.
Analysing the constraints
Let’s see if the system got our own constraints. We look at n.constraints which combines summarises constraints going into the linear problem
[13]:
n.constraints
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In [13], line 1
----> 1 n.constraints
NameError: name 'n' is not defined
The last three entries show our constraints. As ‘soc_lower_bound’ is time-dependent, the pnl value is set to True.
Let’s check whether out two custom constraint are fulfilled:
[14]:
n.links.loc[["battery_power", "battery_discharge"], ["p_nom_opt"]]
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In [14], line 1
----> 1 n.links.loc[["battery_power", "battery_discharge"], ["p_nom_opt"]]
NameError: name 'n' is not defined
[15]:
n.storage_units_t.state_of_charge
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In [15], line 1
----> 1 n.storage_units_t.state_of_charge
NameError: name 'n' is not defined
[16]:
n.generators_t.p.groupby(n.generators.bus, axis=1).sum().sum() / n.loads_t.p.sum().sum()
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In [16], line 1
----> 1 n.generators_t.p.groupby(n.generators.bus, axis=1).sum().sum() / n.loads_t.p.sum().sum()
NameError: name 'n' is not defined
Looks good! Now, let’s see which dual values were parsed. Therefore we have a look into n.dualvalues
[17]:
n.dualvalues
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In [17], line 1
----> 1 n.dualvalues
NameError: name 'n' is not defined
Again we see the last two entries reflect our constraints (the values in the columns play only a minor role). Having a look what the values are:
[18]:
from pypsa.linopt import get_dual
[19]:
get_dual(n, "StorageUnit", "soc_lower_bound")
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In [19], line 1
----> 1 get_dual(n, "StorageUnit", "soc_lower_bound")
NameError: name 'n' is not defined
[20]:
get_dual(n, "battery_discharge", "fixratio")
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In [20], line 1
----> 1 get_dual(n, "battery_discharge", "fixratio")
NameError: name 'n' is not defined
[21]:
get_dual(n, "Bus", "production_share")
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In [21], line 1
----> 1 get_dual(n, "Bus", "production_share")
NameError: name 'n' is not defined
Side note
Some of the predefined constraints are stored in components itself like n.lines_t.mu_upper or n.buses_t.marginal_price, this is the case if their are designated columns are spots for those. All other dual are under the hook stored in n.duals