Problem with formulating lhs constraint in linopy (flow based market coupling constraint)

[gktho...@gmail.com]

Dear all,

I am currently trying to implement a formulation for flow based market coupling, mainly the constraint which says that the flow of a line is equal to the zonal PTDFs multiplied with the zonal net positions.

The issue lies 100 % with my inability in constructing the data arrays correctly, so that it sums up the flows in and out of a zone (==net position) and multiplies it with the zonal ptdf. Any help is therefore highly appreciate it, on how to construct the data so that the constraint works eventually.

I wrote some example code, which I hope, makes the problem reproducible with a a current version of pypsa-eur:

def normed(s): return s/s.sum()

n = pypsa.Network("../pypsa-eur-playground/pypsa-eur/networks/elec_s_128.nc")
m = n.copy()

internal_lines = n.lines[n.lines.bus0.map(n.buses.country) == n.lines.bus1.map(n.buses.country)].copy()
m.mremove("Line", internal_lines.index)
m.madd("Link", internal_lines.index, bus0=internal_lines.bus0.values, bus1=internal_lines.bus1.values, p_nom=1e5) # constructs a zonal model, in which all internal lines are replaced with links with a very high (non restricting) capacity.

n.determine_network_topology() # these following steps are done to simplify the issue, and look only at one sub network.
reduced_network = n.sub_networks.loc["0", "obj"].buses().index
n = n[reduced_network]
m = m[reduced_network]
n.determine_network_topology()

sn = n.sub_networks.loc["0", "obj"]
sn.calculate_PTDF()

nodal_ptdf = pd.DataFrame(sn.PTDF)
nodal_ptdf.columns = sn.buses_i()
nodal_ptdf.index = sn.lines_i()
nodal_ptdf.columns.name = "zone"

techs = ['biomass', 'onwind', 'ror', 'solar', 'CCGT', 'OCGT', 'coal', 'oil',
         'nuclear', 'offwind-ac', 'offwind-dc', 'lignite',]
zonal_dist = n.generators.query("carrier in @techs").groupby("bus").p_nom.sum()
zone = n.buses.loc[zonal_dist.index, "country"]
gsk = zonal_dist.groupby(zone).transform(lambda s: normed(s))
zonal_ptdf = nodal_ptdf.multiply(gsk).groupby(n.buses.loc[nodal_ptdf.columns, "country"],axis=1).sum()

m.optimize.create_model();
o = m.model
line_s = o["Line-s"]

grouper_lines0 = xr.DataArray(
    m.lines.bus0.map(m.buses.country).values,
    dims=["Line"],
    coords=[m.lines.index],
    name = "zone"
)

grouper_lines1 = xr.DataArray(
    m.lines.bus1.map(m.buses.country).values,
    dims=["Line"],
    coords=[m.lines.index],
    name = "zone"
)

net_positions = line_s.groupby(grouper_lines0).sum() - line_s.groupby(grouper_lines1).sum()

zonal_ptdf = xr.DataArray(zonal_ptdf, dims=["Line", "zone"])
lhs = line_s - net_positions*zonal_ptdf # this part is the problematic one
cons = lhs == 0 # I wrote this constraint only to convey what in the end is supposed to be achieved.

Thanks a lot in advance!

Best,

Georg