Note
You can download this example as a Jupyter notebook or start it in interactive mode.
Optimization with Linopy
In PyPSA v0.22, an additional optimization module was instroduced to the package. It is built on Linopy and aims at
performance as we know it from the native PyPSA optimization (
lopfwithpyomo=False)flexibility as we know from the Pyomo implementation
usability as know from pandas/xarray
Linopy is a stand-alone package and works similar to Pyomo but without the memory overhead and much faster. In the long-term we are planning to slowly migrate towards the Linopy-based optimization only. In order to facilitate the transission from the native PyPSA optimization (lopf with pyomo=False), the module pypsa.optimization.compat provides functions similar to pypsa.linopt. Have a look of our migration guide (next notebook).
If you don’t have any code to migrate, we recommend to directly use the linopy functions instead.
For additional information on the Linopy package, have a look at the documentation.
Let’s get started
Now, we demonstrate the behaviour of the optimization with linopy. The core functions for the optimization can be called via the Network.optimize accessor. The accessor is used for creating, solving, modifying the optimization problem. Further it supports to run different optimzation formulations and provides helper functions.
At first, we run the ordinary linearized optimal power flow (LOPF). We then extend the formulation by some additional constraints.
[1]:
import pypsa
import pandas as pd
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
Cell In [1], line 1
----> 1 import pypsa
2 import pandas as pd
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
In order to make the network a bit more interesting, we modify its data: 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 limits,
[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
…and 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
Ordinary Optimization
Per default the optimization based on linopy mimics the well-known n.lopf optimization. We run it by calling the optimize accessor.
[7]:
n.optimize()
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In [7], line 1
----> 1 n.optimize()
NameError: name 'n' is not defined
Compared to the native optimization, we now have a model instance attached to our network. It is a container of all variables, constraints and the objective function. You can modify this as much as you please, by directly adding or deleting variables or constraints etc.
[8]:
n.model
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In [8], line 1
----> 1 n.model
NameError: name 'n' is not defined
Modify model, optimize and feed back to network
When you have a fresh network and you just want to create the model instance, run
[9]:
n.optimize.create_model();
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In [9], line 1
----> 1 n.optimize.create_model();
NameError: name 'n' is not defined
Through the model instance we gain a lot of flexibility. Let’s say for example we want to remove the Kirchhoff Voltage Law constraint, thus convert the model to a transport model. This can be done via
[10]:
n.model.constraints.remove("Kirchhoff-Voltage-Law")
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In [10], line 1
----> 1 n.model.constraints.remove("Kirchhoff-Voltage-Law")
NameError: name 'n' is not defined
Now, we want to optimize the altered model and feed to solution back to the network. Here again, we use the optimize accessor.
[11]:
n.optimize.solve_model()
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In [11], line 1
----> 1 n.optimize.solve_model()
NameError: name 'n' is not defined
Here we followed the recommended way to run altered models:
Create the model instance -
n.optimize.create_model()Modify the model to your needs
Solve and feed back -
n.optimize.solve_model()
For compatibility reasons the optimize function, also allows to pass a extra_funcionality argument, as we know it from the lopf function. The above behaviour with use of the extra functionality is obtained through
[12]:
def remove_kvl(n, sns):
print("KVL removed!")
n.model.constraints.remove("Kirchhoff-Voltage-Law")
n.optimize(extra_functionality=remove_kvl)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In [12], line 6
2 print("KVL removed!")
3 n.model.constraints.remove("Kirchhoff-Voltage-Law")
----> 6 n.optimize(extra_functionality=remove_kvl)
NameError: name 'n' is not defined
Additional constraints
In the following we examplarily present a set of additional constraints. Note, the dual values of the additional constraints won’t be stored in default data fields in the PyPSA network. But in any case they are stored in the linopy.Model.
Again, we first build the optimization model, add our constraints and finally solve the network. For the first step we use again our accessor optimize to access the function create_model. This returns the linopy model that we can modify.
[13]:
m = n.optimize.create_model() # the return value is the model, let's use it directly!
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In [13], line 1
----> 1 m = n.optimize.create_model() # the return value is the model, let's use it directly!
NameError: name 'n' is not defined
Minimum for state of charge
Assume we want to set a minimum state of charge of 50 MWh in our storage unit. This is done by:
[14]:
sus = m.variables["StorageUnit-state_of_charge"]
m.add_constraints(sus >= 50, name="StorageUnit-minimum_soc")
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In [14], line 1
----> 1 sus = m.variables["StorageUnit-state_of_charge"]
2 m.add_constraints(sus >= 50, name="StorageUnit-minimum_soc")
NameError: name 'm' is not defined
The return value of the add_constraints function is a array with the labels of the constraints. You can access the constraint now through:
[15]:
m.constraints["StorageUnit-minimum_soc"]
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In [15], line 1
----> 1 m.constraints["StorageUnit-minimum_soc"]
NameError: name 'm' is not defined
and inspects its attributes like lhs, sign and rhs, e.g.
[16]:
m.constraints["StorageUnit-minimum_soc"].rhs
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In [16], line 1
----> 1 m.constraints["StorageUnit-minimum_soc"].rhs
NameError: name 'm' is not defined
Fix the ratio between ingoing and outgoing capacity of the Store
The battery in our system is modelled with two links and a store. We should make sure that its charging and discharging capacities, meaning their links, are somehow coupled.
[17]:
capacity = m.variables["Link-p_nom"]
eff = n.links.at["battery_power", "efficiency"]
lhs = capacity["battery_power"] - eff * capacity["battery_discharge"]
m.add_constraints(lhs == 0, name="Link-battery_fix_ratio")
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In [17], line 1
----> 1 capacity = m.variables["Link-p_nom"]
2 eff = n.links.at["battery_power", "efficiency"]
3 lhs = capacity["battery_power"] - eff * capacity["battery_discharge"]
NameError: name 'm' is not defined
Every bus must in total produce the 20% of the total demand
For this, we use the linopy function groupby_sum which follows the pattern from pandas/xarray’s groupby function.
[18]:
total_demand = n.loads_t.p_set.sum().sum()
prod_per_bus = m.variables["Generator-p"].groupby_sum(n.generators.bus).sum("snapshot")
m.add_constraints(prod_per_bus >= total_demand / 5, name="Bus-minimum_production_share")
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In [18], line 1
----> 1 total_demand = n.loads_t.p_set.sum().sum()
2 prod_per_bus = m.variables["Generator-p"].groupby_sum(n.generators.bus).sum("snapshot")
3 m.add_constraints(prod_per_bus >= total_demand / 5, name="Bus-minimum_production_share")
NameError: name 'n' is not defined
[19]:
con = prod_per_bus >= total_demand / 5
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In [19], line 1
----> 1 con = prod_per_bus >= total_demand / 5
NameError: name 'prod_per_bus' is not defined
[20]:
con.lhs
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In [20], line 1
----> 1 con.lhs
NameError: name 'con' is not defined
… and now let’s solve the network again.
[21]:
n.optimize.solve_model()
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In [21], line 1
----> 1 n.optimize.solve_model()
NameError: name 'n' is not defined
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
[22]:
n.model.constraints
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In [22], line 1
----> 1 n.model.constraints
NameError: name 'n' is not defined
The last three entries show our constraints. Let’s check whether out two custom constraint are fulfilled:
[23]:
n.links.loc[["battery_power", "battery_discharge"], ["p_nom_opt"]]
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In [23], line 1
----> 1 n.links.loc[["battery_power", "battery_discharge"], ["p_nom_opt"]]
NameError: name 'n' is not defined
[24]:
n.storage_units_t.state_of_charge
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In [24], line 1
----> 1 n.storage_units_t.state_of_charge
NameError: name 'n' is not defined
[25]:
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 [25], 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.model.dual
[26]:
n.model.dual
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In [26], line 1
----> 1 n.model.dual
NameError: name 'n' is not defined
[27]:
n.model.dual["StorageUnit-minimum_soc"]
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In [27], line 1
----> 1 n.model.dual["StorageUnit-minimum_soc"]
NameError: name 'n' is not defined
[28]:
n.model.dual["Link-battery_fix_ratio"]
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In [28], line 1
----> 1 n.model.dual["Link-battery_fix_ratio"]
NameError: name 'n' is not defined
[29]:
n.model.dual["Bus-minimum_production_share"]
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In [29], line 1
----> 1 n.model.dual["Bus-minimum_production_share"]
NameError: name 'n' is not defined
These are the basic functionalities of the optimize accessor. There are many more functions like abstract optimziation formulations (security constraint optimization, iterative transmission expansion optimization, etc.) or helper functions (fixing optimized capacities, adding load shedding). Try them out if you want!
[30]:
print("\n".join([func for func in n.optimize.__dir__() if not func.startswith("_")]))
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In [30], line 1
----> 1 print("\n".join([func for func in n.optimize.__dir__() if not func.startswith("_")]))
NameError: name 'n' is not defined