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

[2]:

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

WARNING:pypsa.io:Importing network from PyPSA version v0.17.1 while current version is v0.22.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


Modify the network a bit: We set gas generators to non-extendable

[3]:

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


[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


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

"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


[6]:

n.add("Bus", "storebus", carrier="hydro", x=-5, y=55)
["battery_power", "battery_discharge"],
"",
bus0=["Manchester", "storebus"],
bus1=["storebus", "Manchester"],
p_nom=100,
efficiency=0.9,
p_nom_extendable=True,
p_nom_max=1000,
)
"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,
);


## 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.

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

1. Fix the ratio between ingoing and outgoing capacity of the Store

[9]:

def fix_link_cap_ratio(n, snapshots):
eff = n.links.at["battery_power", "efficiency"]
lhs = linexpr(
)
define_constraints(n, lhs, "=", 0, "battery_discharge", attr="fixratio")

1. 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):
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_bus_production(n, snapshots)


…and run the lopf with pyomo=False

[12]:

n.lopf(
pyomo=False,
extra_functionality=extra_functionalities,
)

INFO:pypsa.linopf:Prepare linear problem
INFO:pypsa.linopf:Total preparation time: 0.3s
INFO:pypsa.linopf:Solve linear problem using Glpk solver
INFO:pypsa.linopf:Optimization successful. Objective value: 1.43e+07

[12]:

('ok', 'optimal')


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

[13]:

pnl specification
component name
Generator mu_lower True non_ext, p
mu_upper True non_ext, p
Line mu_upper True s
mu_lower True s
Link mu_upper True p
mu_lower True p
Store mu_upper True e
mu_lower True e
StorageUnit mu_upper True p_dispatch, p_store, state_of_charge
mu_lower True p_dispatch, p_store, state_of_charge
mu_state_of_charge_set True
Generator mu_ramp_limit_up True nonext.
mu_ramp_limit_down True nonext.
StorageUnit mu_state_of_charge True
Store mu_state_of_charge True
SubNetwork mu_kirchhoff_voltage_law True
Bus marginal_price True
GlobalConstraint mu False co2_limit
StorageUnit soc_lower_bound True
battery_discharge fixratio False
Bus production_share False

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

[14]:

p_nom_opt
battery_power 900.0
battery_discharge 1000.0
[15]:

n.storage_units_t.state_of_charge

[15]:

StorageUnit su su2
snapshot
2015-01-01 00:00:00 1835.74 1000.000
2015-01-01 01:00:00 1326.16 1000.000
2015-01-01 02:00:00 1376.16 1000.000
2015-01-01 03:00:00 1426.16 1000.000
2015-01-01 04:00:00 1986.06 1000.000
2015-01-01 05:00:00 50.00 50.000
2015-01-01 06:00:00 50.00 50.000
2015-01-01 07:00:00 100.00 156.727
2015-01-01 08:00:00 50.00 140.120
2015-01-01 09:00:00 50.00 50.000
[16]:

n.generators_t.p.groupby(n.generators.bus, axis=1).sum().sum() / n.loads_t.p.sum().sum()

[16]:

bus
Frankfurt     0.200000
Manchester    0.200000
Norway        0.637047
dtype: float64


Looks good! Now, let’s see which dual values were parsed. Therefore we have a look into n.dualvalues

[17]:

n.dualvalues

[17]:

in_comp pnl
component name
StorageUnit mu_upper True True
mu_lower True True
mu_state_of_charge_set True True
mu_state_of_charge True True
Bus marginal_price True True
StorageUnit soc_lower_bound True True
battery_discharge fixratio False False
Bus production_share True False

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

[19]:

StorageUnit su su2
snapshot
2015-01-01 00:00:00 -0.00000 -0.00000
2015-01-01 01:00:00 -0.00000 -0.00000
2015-01-01 02:00:00 -0.00000 -0.00000
2015-01-01 03:00:00 -0.00000 -0.00000
2015-01-01 04:00:00 -0.00000 -0.00000
2015-01-01 05:00:00 -0.00000 -0.00000
2015-01-01 06:00:00 -157.42300 -4.44444
2015-01-01 07:00:00 -0.00000 -0.00000
2015-01-01 08:00:00 -0.00000 -0.00000
2015-01-01 09:00:00 -5.55556 -67.84860
[20]:

get_dual(n, "battery_discharge", "fixratio")

[20]:

0    567.89
Name: fixratio, dtype: float64

[21]:

get_dual(n, "Bus", "production_share")

[21]:

Bus
London            NaN
Norwich           NaN
Norwich DC        NaN
Manchester   -45.4185
Bremen            NaN
Bremen DC         NaN
Frankfurt    -35.2007
Norway        -0.0000
Norway DC         NaN
storebus          NaN
Name: production_share, dtype: float64


### 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