Example on defining "New Components" compatible with the new Optimization class using Linopy

Question

In old versions of PyPSA (before v0.25), there was an example on defining new components, that was a combined heat and power system (link). However, by introduction of the new optimization tool using Linopy, this example became outdated and then removed after v0.25. My first question is whether or not there has been any new example for defining new components compatible with the new network.optimize() command instead of the old network.lopf().

I tried to define my own component in a similar way by adding two links along with the constraints in the optimize() method under the Network class (instead of lopf()). However, overriding optimize() which belongs to OptimizeAccessor class was not successful. I then tried to add a new method for the Network class (see below). This method should be called just before calling network.optimize(). This method somehow does the work, but I'm not sure if it's a good solution. I wonder if you could suggest any improvement or an alternative way to add a new component (such as this CHP).

import numpy as np
import pandas as pd
import pypsa
from pypsa.descriptors import Dict

# Copy existing components and add new components [Directly copied from the old example]
override_components = pypsa.components.components.copy()

override_component_attrs = Dict(
    {k: v.copy() for k, v in pypsa.components.component_attrs.items()}
)

## Adding CHP new element [Directly copied from the old example]
override_components.loc["CHP"] = ["chps", "Combined heat and power plant.", np.nan]

override_component_attrs["CHP"] = pd.DataFrame(
    columns=["type", "unit", "default", "description", "status"]
)

override_component_attrs["CHP"].loc["name"] = [
    "string",
    "n/a",
    "n/a",
    "Unique name",
    "Input (required)",
]
override_component_attrs["CHP"].loc["bus_fuel"] = [
    "string",
    "n/a",
    "n/a",
    "Name of bus where fuel source is.",
    "Input (required)",
]
override_component_attrs["CHP"].loc["bus_elec"] = [
    "string",
    "n/a",
    "n/a",
    "Name of bus where electricity is supplied.",
    "Input (required)",
]
override_component_attrs["CHP"].loc["bus_heat"] = [
    "string",
    "n/a",
    "n/a",
    "Name of bus where heat is supplied.",
    "Input (required)",
]
override_component_attrs["CHP"].loc["p_nom_extendable"] = [
    "boolean",
    "n/a",
    False,
    "",
    "Input (optional)",
]
override_component_attrs["CHP"].loc["marginal_cost"] = [
    "static or series",
    "currency/MWh",
    0.0,
    "Marginal cost of production of 1 MWh.",
    "Input (optional)",
]
override_component_attrs["CHP"].loc["capital_cost"] = [
    "float",
    "EUR/MW",
    0.0,
    "Capital cost per rating of electricity output.",
    "Input (optional)",
]
override_component_attrs["CHP"].loc["eta_elec"] = [
    "float",
    "n/a",
    1.0,
    "Electrical efficiency with no heat output, i.e. in condensing mode",
    "Input (optional)",
]
override_component_attrs["CHP"].loc["c_v"] = [
    "float",
    "n/a",
    1.0,
    "Loss of fuel for each addition of heat",
    "Input (optional)",
]
override_component_attrs["CHP"].loc["c_m"] = [
    "float",
    "n/a",
    1.0,
    "Backpressure ratio",
    "Input (optional)",
]
override_component_attrs["CHP"].loc["p_nom_ratio"] = [
    "float",
    "n/a",
    1.0,
    "Ratio of max heat output to max electrical output; max heat of 500 MWth and max electricity of 1000 MWth means p_nom_ratio is 0.5",
    "Input (optional)",
]


# Defining a new method called "add_elements" under the "Network" class
class Network(pypsa.Network):
    def __init__(self, *args, **kwargs):
        kwargs["override_components"] = override_components
        kwargs["override_component_attrs"] = override_component_attrs
        super().__init__(*args, **kwargs)

    # def a new method, instead of overriding optimize(self, *args, **kwargs):
    def add_elements(self, *args, **kwargs):
        # at this point check that all the extra links are in place for the CHPs
        if not self.chps.empty:
            # for the electricity generation
            self.madd(
                "Link",
                self.chps.index + " electric",
                bus0=self.chps.bus_fuel.values,
                bus1=self.chps.bus_elec.values,
                p_nom_extendable=self.chps.p_nom_extendable.values,
                capital_cost=self.chps.capital_cost.values * self.chps.eta_elec.values,
                efficiency=self.chps.eta_elec.values,
            )

            # for heat generation
            self.madd(
                "Link",
                self.chps.index + " heat",
                bus0=self.chps.bus_fuel.values,
                bus1=self.chps.bus_heat.values,
                p_nom_extendable=self.chps.p_nom_extendable.values,
                efficiency=self.chps.eta_elec.values / self.chps.c_v.values,
            )
            print("_________________ ADDED CHP Links _______________________")

        
        if "extra_functionality" in kwargs:
            kwargs["extra_functionality"] = kwargs.pop("extra_functionality")

        # Consistency check
        self.consistency_check()

        # Creating the model
        m = self.optimize.create_model()

        # Adding constraintse using the Linopy format (not based on the compat class)
        for ii in self.chps.index:
            p_elec = m.variables["Link-p"].loc[:, ii + " electric"]
            p_heat = m.variables["Link-p"].loc[:, ii + " heat"]
            p_nom_elec = m.variables["Link-p_nom"].loc[ii + " electric"]
            p_nom_heat = m.variables["Link-p_nom"].loc[ii + " heat"]
            eta_elec = self.chps.at[ii, "eta_elec"]
            p_nom_ratio = self.chps.at[ii, "p_nom_ratio"]
            c_v = self.chps.at[ii, "c_v"]
            c_m = self.chps.at[ii, "c_m"]

            # 1) ratio of heat and electric production
            constraint1 = (
                eta_elec * p_nom_ratio * p_nom_elec == eta_elec / c_v * p_nom_heat
            )
            m.add_constraints(constraint1, name="Chp-nom_power_" + ii)

            # 2) back pressure
            constraint2 = c_m * eta_elec / c_v * p_heat <= eta_elec * p_elec
            m.add_constraints(constraint2, name="Chp-back_pressure_" + ii)

            # 3) top iso fuel line
            constraint3 = p_heat + p_elec <= p_nom_elec
            m.add_constraints(constraint3, name="Link-isofuelline_" + ii)
            print("_________________SJ: ADDED all three Constraints _______________________")

        print("_________________SJ: CHP Addition COMPLETED _______________________")