Repository: PyPSA/PyPSA

Issue:
Modular investment capacity
<!-- Please do not post usage questions here. Ask them on the PyPSA mailing list: https://groups.google.com/forum/#!forum/pypsa -->

## Describe the feature you'd like to see

It would be interesting to add the modularity in the investment capacity: rather than having a pure LP formulation, enabling the possibility to install a "number" of components.

I reflected on this topic and I see two options, detailed below:
1. (my favourite solution) add a new integer variable n_nom that describes the number of units installed and add a constraint "p_nom (variable) == p_nom (parameter)  * n_nom (variable)", where the latter p_nom shall be the optional parameter set during the "add/madd" function.
2. revise the formulation of the algoritm from p_nom to n_nom. This may be slightly more efficient from a mathematical point of view, but a breaking change. So I wouldn't recommend this option

This feature may be interesting to have and especially option 1 may be easy to implement.
The needed elements may be:
1. add a new `define_integer` function to add integer-valued function, accounting for the LP format to add integer variables, hence using "General" section https://www.ibm.com/docs/en/icos/12.8.0.0?topic=representation-mip-features-in-lp-file-format but also adding a constraint to specify that they are non-negative
2. Add an appropriate constraint to set: "p_nom (variable) == p_nom (parameter)  * n_nom (variable)".
   To avoid adding additional parameters and since the parameter p_nom generally loses value when using extendable, we may use that same parameter to specify the size of each base unit

This simple implementation is consistent with the current PyPSA structure. A follow-up of this could be a revision of the unit commitment to address multi-units of the same type.
This can be done by extending to integer the variable status and simple adjustments.

With @enricogiglio, we could work on this, but comments are welcome to steer the work in a direction that is acceptable.

To ease the review process, the first PR shall not address the revision of the unit commitment with multi-units

Discussion:

Patch:
diff --git a/doc/release_notes.rst b/doc/release_notes.rst
index 637ef29c..b90416ff 100644
--- a/doc/release_notes.rst
+++ b/doc/release_notes.rst
@@ -69,6 +69,12 @@ PyPSA 0.25.1 (27th July 2023)
   those that already have a designated place in the network are assigned.
   (https://github.com/PyPSA/PyPSA/pull/635)
 
+* Added features ``p_nom_mod``, ``e_nom_mod``, and ``s_nom_mod`` to consider
+  capacity modularity of components. When such feature is non-zero and
+  ``exentability`` is true, the component's optimal nominal capacity becomes
+  a multiple of the capacity of the single component. The optimal number of
+  components is ``n_mod`` (such that p_nom_mod*n_mod=p_nom_opt).
+
 **Changes**
 
 * The function ``get_buses_linemap_and_lines`` was deprecated, in favor of
diff --git a/pypsa/component_attrs/generators.csv b/pypsa/component_attrs/generators.csv
index c7316dcb..cfbecfeb 100644
--- a/pypsa/component_attrs/generators.csv
+++ b/pypsa/component_attrs/generators.csv
@@ -4,6 +4,7 @@ bus,string,n/a,n/a,name of bus to which generator is attached,Input (required)
 control,string,n/a,PQ,"P,Q,V control strategy for PF, must be ""PQ"", ""PV"" or ""Slack"".",Input (optional)
 type,string,n/a,n/a,Placeholder for generator type. Not yet implemented.,Input (optional)
 p_nom,float,MW,0,Nominal power for limits in optimization.,Input (optional)
+p_nom_mod,float,MW,0,Nominal power of the generator module.,Input (optional)
 p_nom_extendable,boolean,n/a,False,Switch to allow capacity p_nom to be extended in optimization.,Input (optional)
 p_nom_min,float,MW,0,"If p_nom is extendable in optimization, set its minimum value.",Input (optional)
 p_nom_max,float,MW,inf,"If p_nom is extendable in optimization, set its maximum value (e.g. limited by technical potential).",Input (optional)
@@ -35,6 +36,7 @@ weight,float,n/a,1,"Weighting of a generator. Only used for network clustering."
 p,series,MW,0.,active power at bus (positive if net generation),Output
 q,series,MVar,0.,reactive power (positive if net generation),Output
 p_nom_opt,float,MW,0.,Optimised nominal power.,Output
+n_mod,int,n/a,0,Optimised number of modular generators if ``p_nom_mod`` is non-zero.,Output
 status,series,n/a,1,"Status (1 is on, 0 is off). Only outputted if committable is True.",Output
 mu_upper,series,currency/MWh,0.,Shadow price of upper p_nom limit,Output
 mu_lower,series,currency/MWh,0.,Shadow price of lower p_nom limit,Output
diff --git a/pypsa/component_attrs/lines.csv b/pypsa/component_attrs/lines.csv
index 881f017c..d33f3c47 100644
--- a/pypsa/component_attrs/lines.csv
+++ b/pypsa/component_attrs/lines.csv
@@ -8,6 +8,7 @@ r,float,Ohm,0,"Series resistance, must be non-zero for DC branch in linear power
 g,float,Siemens,0,Shunt conductivity. Shunt admittance is :math:`y = g + jb`.,Input (optional)
 b,float,Siemens,0,Shunt susceptance. If the line has shunt capacitance :math:`C` in Farads then :math:`b = 2\pi f C` where :math:`f` is the frequency in Hertz. Shunt admittance is :math:`y = g + jb`. Ignored if type defined.,Input (optional)
 s_nom,float,MVA,0,Limit of apparent power which can pass through branch.,Input (optional)
+s_nom_mod,float,MWA,0,Limit of apparent power of the line module.,Input (optional)
 s_nom_extendable,boolean,n/a,False,Switch to allow capacity s_nom to be extended in OPF.,Input (optional)
 s_nom_min,float,MVA,0,"If s_nom is extendable in OPF, set its minimum value.",Input (optional)
 s_nom_max,float,MVA,inf,"If s_nom is extendable in OPF, set its maximum value (e.g. limited by potential).",Input (optional)
@@ -33,5 +34,6 @@ b_pu,float,per unit,0,Per unit shunt susceptance calculated by PyPSA from b and
 x_pu_eff,float,per unit,0,"Effective per unit series reactance for linear power flow, calculated by PyPSA from x, tap_ratio for transformers and bus.v_nom.",Output
 r_pu_eff,float,per unit,0,"Effective per unit series resistance for linear power flow, calculated by PyPSA from x, tap_ratio for transformers and bus.v_nom.",Output
 s_nom_opt,float,MVA,0,Optimised capacity for apparent power.,Output
+n_mod,int,n/a,0,Optimised number of modular lines if ``s_nom_mod`` is non-zero.,Output
 mu_lower,series,currency/MVA,0.,Shadow price of lower s_nom limit  -F \leq f. Always non-negative.,Output
 mu_upper,series,currency/MVA,0.,Shadow price of upper s_nom limit f \leq F. Always non-negative.,Output
diff --git a/pypsa/component_attrs/links.csv b/pypsa/component_attrs/links.csv
index 6890f0a7..c8f4f326 100644
--- a/pypsa/component_attrs/links.csv
+++ b/pypsa/component_attrs/links.csv
@@ -8,6 +8,7 @@ efficiency,static or series,per unit,1,Efficiency of power transfer from bus0 to
 build_year,int,year,0,build year,Input (optional)
 lifetime,float,years,inf,lifetime,Input (optional)
 p_nom,float,MW,0,Limit of active power which can pass through link.,Input (optional)
+p_nom_mod,float,MW,0,Limit of active power of the link module.,Input (optional)
 p_nom_extendable,boolean,n/a,False,Switch to allow capacity p_nom to be extended.,Input (optional)
 p_nom_min,float,MW,0,"If p_nom is extendable, set its minimum value.",Input (optional)
 p_nom_max,float,MW,inf,"If p_nom is extendable, set its maximum value (e.g. limited by potential).",Input (optional)
@@ -34,6 +35,7 @@ ramp_limit_shut_down,float,per unit,1,"Maximum decrease at shut down, per unit o
 p0,series,MW,0.,Active power at bus0 (positive if branch is withdrawing power from bus0).,Output
 p1,series,MW,0.,Active power at bus1 (positive if branch is withdrawing power from bus1).,Output
 p_nom_opt,float,MW,0,Optimised capacity for active power.,Output
+n_mod,int,n/a,0,Optimised number of modular links if ``p_nom_mod`` is non-zero.,Output
 status,series,n/a,1.,"Status (1 is on, 0 is off). Only outputted if committable is True.",Output
 mu_lower,series,currency/MW,0.,Shadow price of lower p_nom limit  -F \leq f. Always non-negative.,Output
 mu_upper,series,currency/MW,0.,Shadow price of upper p_nom limit f \leq F. Always non-negative.,Output
diff --git a/pypsa/component_attrs/storage_units.csv b/pypsa/component_attrs/storage_units.csv
index a053f313..17cc6c54 100644
--- a/pypsa/component_attrs/storage_units.csv
+++ b/pypsa/component_attrs/storage_units.csv
@@ -4,6 +4,7 @@ bus,string,n/a,n/a,Name of bus to which storage unit is attached.,Input (require
 control,string,n/a,PQ,"P,Q,V control strategy for PF, must be ""PQ"", ""PV"" or ""Slack"".",Input (optional)
 type,string,n/a,n/a,Placeholder for storage unit type. Not yet implemented.,Input (optional)
 p_nom,float,MW,0,Nominal power for limits in OPF.,Input (optional)
+p_nom_mod,float,MW,0,Nominal power of the storage unit module.,Input (optional)
 p_nom_extendable,boolean,n/a,False,Switch to allow capacity p_nom to be extended in OPF.,Input (optional)
 p_nom_min,float,MW,0,"If p_nom is extendable in OPF, set its minimum value.",Input (optional)
 p_nom_max,float,MW,inf,"If p_nom is extendable in OPF, set its maximum value (e.g. limited by potential).",Input (optional)
@@ -35,6 +36,7 @@ q,series,MVar,0,reactive power (positive if net generation),Output
 state_of_charge,series,MWh,NaN,State of charge as calculated by the OPF.,Output
 spill,series,MW,0,Spillage for each snapshot.,Output
 p_nom_opt,float,MW,0,Optimised nominal power.,Output
+n_mod,int,n/a,0,Optimised number of moudaler storage units if ``p_nom_mod`` is non-zero.,Output
 mu_upper,series,currency/MWh,0,Shadow price of upper p_nom limit,Output
 mu_lower,series,currency/MWh,0,Shadow price of lower p_nom limit,Output
 mu_state_of_charge_set,series,currency/MWh,0,Shadow price of fixed state of charge state_of_charge_set,Output
diff --git a/pypsa/component_attrs/stores.csv b/pypsa/component_attrs/stores.csv
index daaf97ef..252d2826 100644
--- a/pypsa/component_attrs/stores.csv
+++ b/pypsa/component_attrs/stores.csv
@@ -4,6 +4,7 @@ bus,string,n/a,n/a,Name of bus to which store is attached.,Input (required)
 type,string,n/a,n/a,Placeholder for store type. Not yet implemented.,Input (optional)
 carrier,string,n/a,"","Energy carrier of the Store, e.g. ""heat"" or ""gas"".","Input (optional)"
 e_nom,float,MWh,0,Nominal energy capacity.,Input (optional)
+e_nom_mod,float,MWh,0,Nominal energy capacity of the store module.,Input (optional)
 e_nom_extendable,boolean,n/a,False,Switch to allow capacity e_nom to be extended in OPF.,Input (optional)
 e_nom_min,float,MWh,0,"If e_nom is extendable in OPF, set its minimum value.",Input (optional)
 e_nom_max,float,MWh,inf,"If e_nom is extendable in OPF, set its maximum value (e.g. limited by technical potential).",Input (optional)
@@ -26,6 +27,7 @@ p,series,MW,0,active power at bus (positive if net generation),Output
 q,series,MVar,0,reactive power (positive if net generation),Output
 e,series,MWh,0,Energy as calculated by the OPF.,Output
 e_nom_opt,float,MWh,0,Optimised nominal energy capacity outputed by OPF.,Output
+n_mod,int,n/a,0,Optimised number of modular stores if ``e_nom_mod`` is non-zero.,Output
 mu_upper,series,currency/MWh,0,Shadow price of upper e_nom limit,Output
 mu_lower,series,currency/MWh,0,Shadow price of lower e_nom limit,Output
 mu_energy_balance,series,currency/MWh,0,Shadow price of storage consistency equations,Output
diff --git a/pypsa/component_attrs/transformers.csv b/pypsa/component_attrs/transformers.csv
index 06e3e6a7..fa05c521 100644
--- a/pypsa/component_attrs/transformers.csv
+++ b/pypsa/component_attrs/transformers.csv
@@ -9,6 +9,7 @@ r,float,per unit,0,"Series resistance (per unit, using s_nom as base power); mus
 g,float,per unit,0,"Shunt conductivity (per unit, using s_nom as base power). Ignored if type defined.",Input (optional)
 b,float,per unit,0,"Shunt susceptance (per unit, using s_nom as base power). Ignored if type defined.",Input (optional)
 s_nom,float,MVA,0,Limit of apparent power which can pass through branch.,Input (optional)
+s_nom_mod,float,MWA,0,Limit of apparent power of the transformer module.,Input (optional)
 s_nom_extendable,boolean,n/a,False,Switch to allow capacity s_nom to be extended in OPF.,Input (optional)
 s_nom_min,float,MVA,0,"If s_nom is extendable in OPF, set its minimum value.",Input (optional)
 s_nom_max,float,MVA,inf,"If s_nom is extendable in OPF, set its maximum value (e.g. limited by potential).",Input (optional)
@@ -35,5 +36,6 @@ b_pu,float,per unit,0,Per unit shunt susceptance calculated by PyPSA from b and
 x_pu_eff,float,per unit,0,"Effective per unit series reactance for linear power flow, calculated by PyPSA from x, tap_ratio for transformers and bus.v_nom.",Output
 r_pu_eff,float,per unit,0,"Effective per unit series resistance for linear power flow, calculated by PyPSA from x, tap_ratio for transformers and bus.v_nom.",Output
 s_nom_opt,float,MVA,0,Optimised capacity for apparent power.,Output
+n_mod,int,n/a,0,Optimised number of modular transformers if ``s_nom_mod`` is non-zero.,Output
 mu_lower,series,currency/MVA,0,Shadow price of lower s_nom limit  -F \leq f. Always non-negative.,Output
 mu_upper,series,currency/MVA,0,Shadow price of upper s_nom limit f \leq F. Always non-negative.,Output
diff --git a/pypsa/linopf.py b/pypsa/linopf.py
index 8d3af050..d52b19a2 100644
--- a/pypsa/linopf.py
+++ b/pypsa/linopf.py
@@ -35,6 +35,7 @@ from pypsa.linopt import (
     align_with_static_component,
     define_binaries,
     define_constraints,
+    define_integer,
     define_variables,
     get_con,
     get_var,
@@ -257,6 +258,62 @@ def define_unit_commitment_status_variables(n, sns, c):
     define_binaries(n, (sns, com_i), c, "status", mask=active)
 
 
+def define_modular_variables(n, c, attr):
+    """
+    Initializes variables 'attr' for a given component c to allow a modular
+    expansion of the attribute 'attr_nom' It allows to define 'n_opt', the
+    optimal number of installed modules.
+
+    Parameters
+    ----------
+    n : pypsa.Network
+    c : str
+        network component of which the nominal capacity should be defined
+    attr : str
+        name of the variable to be handled attached to modular constraints, e.g. 'p_nom'
+    """
+    mod_i = n.df(c).query(f"{attr}_extendable and ({attr}_mod>0)").index
+    mod_i = mod_i.rename(f"{c}-ext")
+
+    if (mod_i).empty:
+        return
+
+    define_integer(n, 0, np.inf, c, attr + "-n_mod", axes=mod_i)
+
+
+def define_modular_constraints(n, c, attr):
+    """
+    Sets constraints for fixing modular variables of a given component. It
+    allows to define optimal capacity of a component as multiple of the nominal
+    capacity of the single module.
+
+    Parameters
+    ----------
+    n : pypsa.Network
+    c : str
+        name of the network component
+    attr : str
+        name of the parameter representing the capacity of each module, e.g. 'p_nom'
+    """
+    mod_i = n.df(c).query(f"{attr}_extendable and ({attr}_mod>0)").index
+
+    if (mod_i).empty:
+        return
+
+    modularity = get_var(n, c, attr + "-n_mod")[mod_i]
+    modular_capacity = n.df(c)[f"{attr}_mod"].loc[mod_i]
+    capacity = get_var(n, c, attr)[mod_i]
+
+    lhs = linexpr(
+        (
+            modular_capacity,
+            modularity,
+        ),
+        (-1, capacity),
+    )
+    define_constraints(n, lhs, "=", 0, c, f"mu_{attr}_mod")
+
+
 def define_loss_variables(n, sns, c, transmission_losses):
     branches_i = n.df(c).index
     if branches_i.empty or not transmission_losses:
@@ -1183,23 +1240,29 @@ def prepare_lopf(
     fdo, objective_fn = mkstemp(".txt", "pypsa-objectve-", **tmpkwargs)
     fdc, constraints_fn = mkstemp(".txt", "pypsa-constraints-", **tmpkwargs)
     fdb, bounds_fn = mkstemp(".txt", "pypsa-bounds-", **tmpkwargs)
+    fdg, generals_fn = mkstemp(".txt", "pypsa-generals-", **tmpkwargs)
     fdi, binaries_fn = mkstemp(".txt", "pypsa-binaries-", **tmpkwargs)
     fdp, problem_fn = mkstemp(".lp", "pypsa-problem-", **tmpkwargs)
 
     n.objective_f = open(objective_fn, mode="w")
     n.constraints_f = open(constraints_fn, mode="w")
     n.bounds_f = open(bounds_fn, mode="w")
+    n.generals_f = open(generals_fn, mode="w")
     n.binaries_f = open(binaries_fn, mode="w")
 
     n.objective_f.write("\ LOPF \n\nmin\nobj:\n")
     n.constraints_f.write("\n\ns.t.\n\n")
     n.bounds_f.write("\nbounds\n")
+    n.generals_f.write("\ngeneral\n")
     n.binaries_f.write("\nbinary\n")
 
     for c in n.passive_branch_components:
         define_loss_variables(n, snapshots, c, transmission_losses)
     for c, attr in lookup.query("nominal and not handle_separately").index:
         define_nominal_for_extendable_variables(n, c, attr)
+        define_modular_variables(n, c, attr)
+        define_modular_constraints(n, c, attr)
+
         # define constraint for newly installed capacity per investment period
         define_growth_limit(n, snapshots, c, attr)
         # define_fixed_variable_constraints(n, snapshots, c, attr, pnl=False)
@@ -1249,6 +1312,7 @@ def prepare_lopf(
         ("bounds_f", fdb),
         ("constraints_f", fdc),
         ("objective_f", fdo),
+        ("generals_f", fdg),
         ("binaries_f", fdi),
     ):
         getattr(n, f).close()
@@ -1257,7 +1321,7 @@ def prepare_lopf(
 
     # concatenate files
     with open(problem_fn, "wb") as wfd:
-        for f in [objective_fn, constraints_fn, bounds_fn, binaries_fn]:
+        for f in [objective_fn, constraints_fn, bounds_fn, generals_fn, binaries_fn]:
             with open(f, "rb") as fd:
                 shutil.copyfileobj(fd, wfd)
             if not keep_files:
@@ -1329,7 +1393,10 @@ def assign_solution(
                 non_ext = n.df(c)[attr]
                 n.df(c)[attr + "_opt"] = sol.reindex(non_ext.index).fillna(non_ext)
             else:
-                n.sols[c].df[attr] = sol
+                if attr.endswith("-n_mod"):
+                    n.df(c)["n_mod"].update(sol)
+                else:
+                    n.sols[c].df[attr] = sol
 
     n.sols = Dict()
     n.solutions = pd.DataFrame(index=n.variables.index, columns=["in_comp", "pnl"])
diff --git a/pypsa/linopt.py b/pypsa/linopt.py
index 8cf1e5de..7deddcf0 100644
--- a/pypsa/linopt.py
+++ b/pypsa/linopt.py
@@ -138,6 +138,34 @@ def define_binaries(n, axes, name, attr="", spec="", mask=None):
     return var
 
 
+def define_integer(n, lower, upper, name, attr="", axes=None, spec=""):
+    """
+    Defines integer-variable(s) for pypsa-network. The variables are stored in
+    the network object under n.vars with key of the variable name.
+
+    Parameters
+    ----------
+    n : pypsa.Network
+    lower, upper: lower bound (lb) and upper bound (ub) and pass it to define_integer
+    name : str
+        general name of the variable (or component which the variable is
+        referring to). The variable will then be stored under:
+
+            * n.vars[name].pnl if the variable is two-dimensional
+            * n.vars[name].df if the variable is one-dimensional
+    attr : str default ''
+        Specifying name of the variable, defines under which name the variable(s)
+        are stored in n.vars[name].pnl if two-dimensional or in n.vars[name].df
+        if one-dimensional. e.g. 'n_opt'
+    axes : pd.Index or tuple of pd.Index objects
+        Specifies the axes and therefore the shape of the variables.
+    """
+    var = write_integer(n, lower, upper, axes, mask=None)
+
+    set_varref(n, var, name, attr, spec=spec)
+    return var
+
+
 def define_constraints(
     n, lhs, sense, rhs, name, attr="", axes=None, spec="", mask=None
 ):
@@ -292,6 +320,30 @@ def write_binary(n, axes, mask=None):
     return to_pandas(variables, *axes)
 
 
+def write_integer(n, lower, upper, axes, mask=None):
+    """
+    Writer function for writing out integer-variables.
+
+    According to the axes it writes out integer for each entry the
+    pd.DataFrame spanned by axes. Returns a frame with variable
+    references.
+    """
+    axes, shape, size = _get_handlers(axes)
+    n._xCounter += size
+    variables = np.arange(n._xCounter - size, n._xCounter).reshape(shape)
+    exprs = "x" + _str_array(variables, True) + "\n"
+    if mask is not None:
+        exprs = np.where(mask, exprs, "")
+        variables = np.where(mask, variables, -1)
+    n.generals_f.write(join_exprs(exprs))
+    if not size:
+        return pd.Series(dtype=float)
+    lower, upper = _str_array(lower), _str_array(upper)
+    exprs = lower + " <= x" + _str_array(variables, True) + " <= " + upper + "\n"
+    n.bounds_f.write(join_exprs(exprs))
+    return to_pandas(variables, *axes)
+
+
 def write_objective(n, terms):
     """
     Writer function for writing out one or multiple objective terms.
diff --git a/pypsa/optimization/constraints.py b/pypsa/optimization/constraints.py
index 47ee7844..9d7c2775 100644
--- a/pypsa/optimization/constraints.py
+++ b/pypsa/optimization/constraints.py
@@ -661,6 +661,34 @@ def define_fixed_nominal_constraints(n, c, attr):
     n.model.add_constraints(var, "=", fix, f"{c}-{attr}_set")
 
 
+def define_modular_constraints(n, c, attr):
+    """
+    Sets constraints for fixing modular variables of a given component. It
+    allows to define optimal capacity of a component as multiple of the nominal
+    capacity of the single module.
+
+    Parameters
+    ----------
+    n : pypsa.Network
+    c : str
+        name of the network component
+    attr : str
+        name of the variable, e.g. 'n_opt'
+    """
+    m = n.model
+    mod_i = n.df(c).query(f"{attr}_extendable and ({attr}_mod>0)").index
+
+    if (mod_i).empty:
+        return
+
+    modularity = m.variables[f"{c}-n_mod"]
+    modular_capacity = n.df(c)[f"{attr}_mod"].loc[mod_i]
+    capacity = m.variables[f"{c}-{attr}"].loc[mod_i]
+
+    con = capacity - modularity * modular_capacity.values == 0
+    n.model.add_constraints(con, name=f"{c}-{attr}_modularity", mask=None)
+
+
 def define_fixed_operation_constraints(n, sns, c, attr):
     """
     Sets constraints for fixing time-dependent variables of a given component
diff --git a/pypsa/optimization/optimize.py b/pypsa/optimization/optimize.py
index e0207113..1b5cfd63 100644
--- a/pypsa/optimization/optimize.py
+++ b/pypsa/optimization/optimize.py
@@ -27,6 +27,7 @@ from pypsa.optimization.constraints import (
     define_fixed_operation_constraints,
     define_kirchhoff_voltage_constraints,
     define_loss_constraints,
+    define_modular_constraints,
     define_nodal_balance_constraints,
     define_nominal_constraints_for_extendables,
     define_operational_constraints_for_committables,
@@ -47,6 +48,7 @@ from pypsa.optimization.global_constraints import (
 )
 from pypsa.optimization.variables import (
     define_loss_variables,
+    define_modular_variables,
     define_nominal_variables,
     define_operational_variables,
     define_shut_down_variables,
@@ -241,6 +243,7 @@ def create_model(
     # Define variables
     for c, attr in lookup.query("nominal").index:
         define_nominal_variables(n, c, attr)
+        define_modular_variables(n, c, attr)
 
     for c, attr in lookup.query("not nominal and not handle_separately").index:
         define_operational_variables(n, sns, c, attr)
@@ -259,6 +262,7 @@ def create_model(
     for c, attr in lookup.query("nominal").index:
         define_nominal_constraints_for_extendables(n, c, attr)
         define_fixed_nominal_constraints(n, c, attr)
+        define_modular_constraints(n, c, attr)
 
     for c, attr in lookup.query("not nominal and not handle_separately").index:
         define_operational_constraints_for_non_extendables(
@@ -349,6 +353,8 @@ def assign_solution(n):
 
             else:
                 set_from_frame(n, c, attr, df)
+        elif attr == "n_mod":
+            n.df(c)[attr].update(df)
         else:
             n.df(c)[attr + "_opt"].update(df)
 
diff --git a/pypsa/optimization/variables.py b/pypsa/optimization/variables.py
index 2770ee94..6713bd76 100644
--- a/pypsa/optimization/variables.py
+++ b/pypsa/optimization/variables.py
@@ -97,6 +97,29 @@ def define_nominal_variables(n, c, attr):
     n.model.add_variables(coords=[ext_i], name=f"{c}-{attr}")
 
 
+def define_modular_variables(n, c, attr):
+    """
+    Initializes variables 'attr' for a given component c to allow a modular
+    expansion of the attribute 'attr_nom' It allows to define 'n_opt', the
+    optimal number of installed modules.
+
+    Parameters
+    ----------
+    n : pypsa.Network
+    c : str
+        network component of which the nominal capacity should be defined
+    attr : str
+        name of the variable to be handled attached to modular constraints, e.g. 'p_nom'
+    """
+    mod_i = n.df(c).query(f"{attr}_extendable and ({attr}_mod>0)").index
+    mod_i = mod_i.rename(f"{c}-ext")
+
+    if (mod_i).empty:
+        return
+
+    n.model.add_variables(lower=0, coords=[mod_i], name=f"{c}-n_mod", integer=True)
+
+
 def define_spillage_variables(n, sns):
     """
     Defines the spillage variables for storage units.
diff --git a/pypsa/variables.csv b/pypsa/variables.csv
index d0dadae6..a4951889 100644
--- a/pypsa/variables.csv
+++ b/pypsa/variables.csv
@@ -4,20 +4,26 @@ Generator,status,False,False,True
 Generator,start_up,False,False,True
 Generator,shut_down,False,False,True
 Generator,p_nom,False,True,False
+Generator,n_mod,False,False,True
 Line,s,False,False,False
 Line,s_nom,False,True,False
+Line,n_mod,False,False,True
 Transformer,s,False,False,False
 Transformer,s_nom,False,True,False
+Transformer,n_mod,False,False,True
 Link,p,True,False,False
 Link,status,False,False,True
 Link,start_up,False,False,True
 Link,shut_down,False,False,True
 Link,p_nom,False,True,False
+Link,n_mod,False,False,True
 Store,e,False,False,False
 Store,p,True,False,True
 Store,e_nom,False,True,False
+Store,n_mod,False,False,True
 StorageUnit,p_dispatch,True,False,False
 StorageUnit,p_store,False,False,False
 StorageUnit,state_of_charge,False,False,False
 StorageUnit,p_nom,False,True,False
+StorageUnit,n_mod,False,False,True
 StorageUnit,spill,False,False,True


Test patch:
diff --git a/test/test_lopf_modularity.py b/test/test_lopf_modularity.py
new file mode 100644
index 00000000..946d1ded
--- /dev/null
+++ b/test/test_lopf_modularity.py
@@ -0,0 +1,55 @@
+# -*- coding: utf-8 -*-
+import numpy as np
+import pandas as pd
+import pytest
+from conftest import SUPPORTED_APIS, optimize
+from numpy.testing import assert_array_almost_equal as equal
+
+import pypsa
+
+MODULARITY_APIS = ["linopy", "native"]
+
+
+@pytest.mark.parametrize("api", MODULARITY_APIS)
+def test_modular_components(api):
+    """
+    This test is based on https://pypsa.readthedocs.io/en/latest/examples/unit-
+    commitment.html and is not very comprehensive.
+    """
+    n = pypsa.Network()
+
+    snapshots = range(4)
+
+    n.set_snapshots(snapshots)
+
+    n.add("Bus", "bus")
+
+    n.add(
+        "Generator",
+        "gas",
+        bus="bus",
+        p_nom_extendable=True,
+        marginal_cost=10,
+        p_nom_max=5000,
+        p_nom_mod=1000,
+    )
+
+    n.add(
+        "Store",
+        "Store_unit",
+        bus="bus",
+        e_nom_extendable=True,
+        e_nom_max=2000,
+        e_nom_mod=100,
+        capital_cost=100,
+    )
+
+    n.add("Load", "load", bus="bus", p_set=[4000, 5000, 6000, 800])
+
+    optimize(n, api)
+
+    expected_n_opt_gen = np.array([5], dtype=float).T
+    expected_n_opt_store = np.array([10], dtype=float).T
+
+    equal(n.generators.n_mod, expected_n_opt_gen)
+    equal(n.stores.n_mod, expected_n_opt_store)