Note

You can download this example as a Jupyter notebook or start it in interactive mode.

Optimization with Linopy - Migrate Extra Functionalities#

The extra funcionalities for the native pypsa optimization code are mostly using the function of the pypsa.linopt model. Here we show how you can recycle large parts of your code by using the compatibility functions from the pypsa.optimization.compat module.

These are

define_variables
define_constraints
get_var
linexpr

You might want to use them if you have extra_functionalities written for the native optimization code. However, expect some hick-ups, as some operations might behave differently.

Let’s import pypsa and the compat functions:

[1]:

import pypsa
from pypsa.optimization.compat import get_var, define_constraints, linexpr

ERROR 1: PROJ: proj_create_from_database: Open of /home/docs/checkouts/readthedocs.org/user_builds/pypsa/conda/latest/share/proj failed

We load the same network from the previous section into memory:

[2]:

n = pypsa.examples.ac_dc_meshed(from_master=True)

n.generators.loc[n.generators.carrier == "gas", "p_nom_extendable"] = False
n.generators.loc[n.generators.carrier == "gas", "ramp_limit_down"] = 0.2
n.generators.loc[n.generators.carrier == "gas", "ramp_limit_up"] = 0.2

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

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,
);

WARNING:pypsa.io:Importing network from PyPSA version v0.17.1 while current version is v0.27.1. Read the release notes at https://pypsa.readthedocs.io/en/latest/release_notes.html to prepare your network for import.

INFO:pypsa.io:Imported network ac-dc-meshed.nc has buses, carriers, generators, global_constraints, lines, links, loads

And define the extra functionalities as we defined them for the native code in here

[3]:

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")

With the compat functions, this will work as expected. Let’s go on to the next one.

[4]:

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.loc["battery_power"]), (-eff, vars_link.loc["battery_discharge"])
    )
    define_constraints(n, lhs, "=", 0, "battery_discharge", attr="fixratio")

This function as well should not make any problems. Let’s go on.

[5]:

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"
    )

Here, we come into difficult terrain. The groupby function won’t work since the linexpr functions returns some sort of a xarray object (a LinearExpression object, derived from xarray.Dataset).

Instead, we have to rewrite parts:

remove the axis argument
use a xaray DataArray as argument
explicitly sum over all snapshots afterwards. This has nothing to do with groupby but with the fact that we want to limit to total production.

[6]:

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.to_xarray())
        .sum()
        .sum()
    )
    define_constraints(
        n, prod_per_bus, ">=", total_demand / 5, "Bus", "production_share"
    )

[7]:

def extra_functionalities(n, snapshots):
    minimal_state_of_charge(n, snapshots)
    fix_link_cap_ratio(n, snapshots)
    fix_bus_production(n, snapshots)

[8]:

n.optimize(
    extra_functionality=extra_functionalities,
)

WARNING:pypsa.components:The following lines have zero x, which could break the linear load flow:
Index(['2', '3', '4'], dtype='object', name='Line')

WARNING:pypsa.components:The following lines have zero r, which could break the linear load flow:
Index(['0', '1', '5', '6'], dtype='object', name='Line')

WARNING:pypsa.components:The following lines have zero x, which could break the linear load flow:
Index(['2', '3', '4'], dtype='object', name='Line')

WARNING:pypsa.components:The following lines have zero r, which could break the linear load flow:
Index(['0', '1', '5', '6'], dtype='object', name='Line')

INFO:linopy.model: Solve problem using Glpk solver

INFO:linopy.io: Writing time: 0.15s

INFO:linopy.solvers:GLPSOL--GLPK LP/MIP Solver 5.0
Parameter(s) specified in the command line:
 --lp /tmp/linopy-problem-k1_ya143.lp --output /tmp/linopy-solve-a78srjv8.sol
Reading problem data from '/tmp/linopy-problem-k1_ya143.lp'...
770 rows, 300 columns, 1653 non-zeros
4350 lines were read
GLPK Simplex Optimizer 5.0
770 rows, 300 columns, 1653 non-zeros
Preprocessing...
546 rows, 288 columns, 1416 non-zeros
Scaling...
 A: min|aij| =  9.693e-03  max|aij| =  2.000e+00  ratio =  2.063e+02
GM: min|aij| =  4.602e-01  max|aij| =  2.173e+00  ratio =  4.721e+00
EQ: min|aij| =  2.125e-01  max|aij| =  1.000e+00  ratio =  4.706e+00
Constructing initial basis...
Size of triangular part is 545
      0: obj =   3.234549052e+03 inf =   8.146e+04 (101)
    145: obj =   2.631753436e+07 inf =   2.842e-14 (0) 1
*   259: obj =   1.415959968e+07 inf =   2.618e-12 (0) 1
OPTIMAL LP SOLUTION FOUND
Time used:   0.0 secs
Memory used: 0.8 Mb (850071 bytes)
Writing basic solution to '/tmp/linopy-solve-a78srjv8.sol'...

INFO:linopy.constants: Optimization successful:
Status: ok
Termination condition: optimal
Solution: 300 primals, 770 duals
Objective: 1.42e+07
Solver model: not available
Solver message: optimal

/home/docs/checkouts/readthedocs.org/user_builds/pypsa/conda/latest/lib/python3.11/site-packages/pypsa/optimization/optimize.py:357: FutureWarning: A value is trying to be set on a copy of a DataFrame or Series through chained assignment using an inplace method.
The behavior will change in pandas 3.0. This inplace method will never work because the intermediate object on which we are setting values always behaves as a copy.

For example, when doing 'df[col].method(value, inplace=True)', try using 'df.method({col: value}, inplace=True)' or df[col] = df[col].method(value) instead, to perform the operation inplace on the original object.


  n.df(c)[attr + "_opt"].update(df)
INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-fix-p-lower, Generator-fix-p-upper, Generator-ext-p-lower, Generator-ext-p-upper, Generator-fix-p-ramp_limit_up, Generator-fix-p-ramp_limit_down, Line-ext-s-lower, Line-ext-s-upper, Link-ext-p-lower, Link-ext-p-upper, Store-ext-e-lower, Store-ext-e-upper, StorageUnit-ext-p_dispatch-lower, StorageUnit-ext-p_dispatch-upper, StorageUnit-ext-p_store-lower, StorageUnit-ext-p_store-upper, StorageUnit-ext-state_of_charge-lower, StorageUnit-ext-state_of_charge-upper, StorageUnit-state_of_charge_set, Kirchhoff-Voltage-Law, StorageUnit-energy_balance, Store-energy_balance, StorageUnit-soc_lower_bound were not assigned to the network.

[8]:

('ok', 'optimal')