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
join_exprs

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, join_exprs, linexpr

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.21.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["battery_power"]), (-eff, vars_link["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:

use groupby_sum function instead groupby
remove the axis 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_sum(n.generators.bus)
        .sum("snapshot")
    )
    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.expressions:Converting group pandas.Series to xarray.DataArray
INFO:linopy.model: Solve linear problem using Glpk solver
INFO:linopy.io: Writing time: 0.41s
INFO:linopy.model: Optimization successful. Objective value: 1.43e+07

GLPSOL--GLPK LP/MIP Solver 5.0
Parameter(s) specified in the command line:
 --lp /tmp/linopy-problem-2p035mi9.lp --output /tmp/linopy-solve-j9usk48p.sol
Reading problem data from '/tmp/linopy-problem-2p035mi9.lp'...
772 rows, 300 columns, 1653 non-zeros
5167 lines were read
GLPK Simplex Optimizer 5.0
772 rows, 300 columns, 1653 non-zeros
Preprocessing...
548 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 547
      0: obj =   3.234550066e+03 inf =   8.754e+04 (102)
    195: obj =   2.782468542e+07 inf =   2.842e-14 (0) 1
*   340: obj =   1.433711360e+07 inf =   4.356e-11 (0) 1
OPTIMAL LP SOLUTION FOUND
Time used:   0.0 secs
Memory used: 0.9 Mb (982675 bytes)
Writing basic solution to '/tmp/linopy-solve-j9usk48p.sol'...

[8]:

('ok', 'optimal')