From a429189bd369569a92174efa71e161b44450be7a Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Pedro=20L=2E=20Magalh=C3=A3es?= <pmlpm@posteo.de>
Date: Sat, 1 Jun 2024 20:05:27 +0200
Subject: [PATCH] Separated network restrictions from the main model module.
 Tests work. Started to work on converters and prices but did not test.

---
 .../problems/esipp/blocks/converters.py       | 1007 +++++++++++++++++
 .../problems/esipp/blocks/networks.py         |  746 ++++++++++++
 src/topupopt/problems/esipp/blocks/prices.py  |  205 ++++
 src/topupopt/problems/esipp/model.py          |  834 +-------------
 tests/test_esipp_network.py                   |   20 +-
 tests/test_esipp_problem.py                   |    6 +-
 6 files changed, 1969 insertions(+), 849 deletions(-)
 create mode 100644 src/topupopt/problems/esipp/blocks/converters.py
 create mode 100644 src/topupopt/problems/esipp/blocks/networks.py
 create mode 100644 src/topupopt/problems/esipp/blocks/prices.py

diff --git a/src/topupopt/problems/esipp/blocks/converters.py b/src/topupopt/problems/esipp/blocks/converters.py
new file mode 100644
index 0000000..ca83d2c
--- /dev/null
+++ b/src/topupopt/problems/esipp/blocks/converters.py
@@ -0,0 +1,1007 @@
+# imports
+
+import pyomo.environ as pyo
+
+# *****************************************************************************
+# *****************************************************************************
+
+def add_converters(
+    model: pyo.AbstractModel,
+    enable_default_values: bool = True,
+    enable_validation: bool = True,
+    enable_initialisation: bool = True,
+):
+
+    # *************************************************************************
+    # *************************************************************************
+
+    # systems
+
+    # set of all systems
+
+    model.set_I = pyo.Set()
+
+    # set of optional systems
+
+    model.set_I_new = pyo.Set(within=model.set_I)
+
+    # *************************************************************************
+
+    # inputs
+
+    # set of inputs (indexed by system)
+
+    model.set_M = pyo.Set(model.set_I)
+
+    # set of inputs modelled using non-negative real variables
+
+    model.set_M_nnr = pyo.Set(model.set_I, within=model.set_M)
+
+    # set of inputs modelled using binary variables
+
+    model.set_M_bin = pyo.Set(model.set_I, within=model.set_M)
+
+    # set of amplitude-constrained inputs
+
+    model.set_M_dim = pyo.Set(model.set_I_new, within=model.set_M)
+
+    # set of amplitude-constrained inputs
+
+    model.set_M_fix = pyo.Set(model.set_I, within=model.set_M)
+
+    # set of externality-inducing inputs
+
+    model.set_M_ext = pyo.Set(model.set_I, within=model.set_M)
+
+    # *************************************************************************
+
+    # outputs
+
+    # set of outputs (indexed by system)
+
+    model.set_R = pyo.Set(model.set_I)
+
+    # set of outputs with fixed bounds
+
+    model.set_R_fix = pyo.Set(model.set_I, within=model.set_R)
+
+    # set of positive amplitude-constrained outputs
+
+    model.set_R_dim_pos = pyo.Set(model.set_I, within=model.set_R)
+
+    # set of negative amplitude-constrained outputs
+
+    model.set_R_dim_neg = pyo.Set(model.set_I, within=model.set_R)
+
+    # set of amplitude-limited outputs with matching pos. and neg. amplitudes
+
+    model.set_R_dim_eq = pyo.Set(model.set_I, within=model.set_R)
+
+    # set of outputs (indexed by system) inducing externalities
+
+    model.set_R_ext = pyo.Set(model.set_I)
+
+    # *************************************************************************
+
+    # states
+
+    # set of states
+
+    model.set_N = pyo.Set(model.set_I)
+
+    # set of states with fixed bounds
+
+    model.set_N_fix = pyo.Set(model.set_I, within=model.set_N)
+
+    # set of positive amplitude-constrained states
+
+    model.set_N_dim_pos = pyo.Set(model.set_I, within=model.set_N)
+
+    # set of negative amplitude-constrained states
+
+    model.set_N_dim_neg = pyo.Set(model.set_I, within=model.set_N)
+
+    # set of amplitude-limited states with matching pos. and neg. amplitudes
+
+    model.set_N_dim_eq = pyo.Set(model.set_I, within=model.set_N)
+
+    # set of states (indexed by system) inducing externalities
+
+    model.set_N_ext = pyo.Set(model.set_I, within=model.set_N)
+
+    # set of positive state variation-penalised states
+
+    model.set_N_pos_var = pyo.Set(model.set_I, within=model.set_N)
+
+    # set of negative state variation-penalised states
+
+    model.set_N_neg_var = pyo.Set(model.set_I, within=model.set_N)
+
+    # set of upper reference violation-penalised states
+
+    model.set_N_ref_u = pyo.Set(model.set_I, within=model.set_N)
+
+    # set of lower reference violation-penalised states
+
+    model.set_N_ref_d = pyo.Set(model.set_I, within=model.set_N)
+
+    # *************************************************************************
+    # *************************************************************************
+
+    # sparse index sets
+
+    # *************************************************************************
+    # *************************************************************************
+
+    # inputs
+
+    # set of IM tuples
+
+    def init_set_IM(m):
+        return ((i, m_i) for i in m.set_I for m_i in m.set_M[i])
+
+    model.set_IM = pyo.Set(dimen=2, initialize=init_set_IM)
+
+    # set of IM tuples for systems with binary signals
+
+    def init_set_IM_bin(m):
+        return ((i, m_i) for (i, m_i) in m.set_IM if m_i in m.set_M_bin[i])
+
+    model.set_IM_bin = pyo.Set(dimen=2, initialize=init_set_IM_bin, within=model.set_IM)
+
+    # set of IM tuples for tech. with dimensionable reference mode levels
+
+    def init_set_IM_dim(m):
+        return ((i, m_i) for (i, m_i) in m.set_IM if m_i in m.set_M_dim[i])
+
+    model.set_IM_dim = pyo.Set(dimen=2, initialize=init_set_IM_dim, within=model.set_IM)
+
+    # set of IM tuples for fixed amplitude inputs
+
+    def init_set_IM_fix(m):
+        return ((i, m_i) for (i, m_i) in m.set_IM if m_i in m.set_M_fix[i])
+
+    model.set_IM_fix = pyo.Set(dimen=2, initialize=init_set_IM_fix, within=model.set_IM)
+
+    # set of IM tuples for technologies whose modes can induce externalities
+
+    def init_set_IM_ext(m):
+        return ((i, m_i) for (i, m_i) in m.set_IM if m_i in m.set_M_ext[i])
+
+    model.set_IM_ext = pyo.Set(dimen=2, initialize=init_set_IM_ext, within=model.set_IM)
+
+    # *************************************************************************
+
+    # states
+
+    # set of IN tuples
+
+    def init_set_IN(m):
+        return (
+            (i, n_i) for i in m.set_I for n_i in m.set_N[i]  # IN tuple
+        )  # for each state
+
+    model.set_IN = pyo.Set(dimen=2, initialize=init_set_IN)
+
+    # set of IN tuples for states with fixed bounds
+
+    def init_set_IN_fix(m):
+        return ((i, n_i) for i in m.set_I for n_i in m.set_N_fix[i])
+
+    model.set_IN_fix = pyo.Set(dimen=2, initialize=init_set_IN_fix)
+
+    # set of IN tuples for converters with amplitude-constrained states
+
+    def init_set_IN_dim_eq(m):
+        return ((i, n_i) for (i, n_i) in m.set_IN if n_i in m.set_N_dim_eq[i])
+
+    model.set_IN_dim_eq = pyo.Set(
+        dimen=2, initialize=init_set_IN_dim_eq, within=model.set_IN
+    )
+
+    # set of IN tuples for converters with pos. amplitude-constrained states
+
+    def init_set_IN_dim_pos(m):
+        return ((i, n_i) for (i, n_i) in m.set_IN if n_i in m.set_N_dim_pos[i])
+
+    model.set_IN_dim_pos = pyo.Set(
+        dimen=2, initialize=init_set_IN_dim_pos, within=model.set_IN
+    )
+
+    # set of IN tuples for converters with neg. amplitude-constrained states
+
+    def init_set_IN_dim_neg(m):
+        return ((i, n_i) for (i, n_i) in m.set_IN if n_i in m.set_N_dim_neg[i])
+
+    model.set_IN_dim_neg = pyo.Set(
+        dimen=2, initialize=init_set_IN_dim_neg, within=model.set_IN
+    )
+
+    # set of IN tuples for converters with externality-inducing states
+
+    def init_set_IN_ext(m):
+        return ((i, n_i) for (i, n_i) in m.set_IN if n_i in m.set_N_ext[i])
+
+    model.set_IN_ext = pyo.Set(dimen=2, initialize=init_set_IN_ext, within=model.set_IN)
+
+    # set of IN tuples for positive variation-penalised states
+
+    def init_set_IN_pos_var(m):
+        return ((i, n_i) for (i, n_i) in m.set_IN if n_i in m.set_N_pos_var[i])
+
+    model.set_IN_pos_var = pyo.Set(
+        dimen=2, initialize=init_set_IN_pos_var, within=model.set_IN
+    )
+
+    # set of IN tuples for negative variation-penalised states
+
+    def init_set_IN_neg_var(m):
+        return ((i, n_i) for (i, n_i) in m.set_IN if n_i in m.set_N_neg_var[i])
+
+    model.set_IN_neg_var = pyo.Set(
+        dimen=2, initialize=init_set_IN_neg_var, within=model.set_IN
+    )
+
+    # set of IN tuples for upper reference violation penalised states
+
+    def init_set_IN_ref_u(m):
+        return ((i, n_i) for (i, n_i) in m.set_IN if n_i in m.set_N_ref_u[i])
+
+    model.set_IN_ref_u = pyo.Set(
+        dimen=2, initialize=init_set_IN_ref_u, within=model.set_IN
+    )
+
+    # set of IN tuples for lower reference violation penalised states
+
+    def init_set_IN_ref_d(m):
+        return ((i, n_i) for (i, n_i) in m.set_IN if n_i in m.set_N_ref_d[i])
+
+    model.set_IN_ref_d = pyo.Set(
+        dimen=2, initialize=init_set_IN_ref_d, within=model.set_IN
+    )
+
+    # *************************************************************************
+
+    # outputs
+
+    # set of IR tuples
+
+    def init_set_IR(m):
+        return ((i, r_i) for i in m.set_I for r_i in m.set_R[i])
+
+    model.set_IR = pyo.Set(dimen=2, initialize=init_set_IR)
+
+    # set of IR tuples for outputs with fixed bounds
+
+    def init_set_IR_fix(m):
+        return ((i, r_i) for i in m.set_I for r_i in m.set_R_fix[i])
+
+    model.set_IR_fix = pyo.Set(dimen=2, initialize=init_set_IR_fix)
+
+    # set of IR tuples for converters with matching pos. and neg. out. amp. limits
+
+    def init_set_IR_dim_eq(m):
+        return ((i, r_i) for (i, r_i) in m.set_IR if r_i in m.set_R_dim_eq[i])
+
+    model.set_IR_dim_eq = pyo.Set(dimen=2, initialize=init_set_IR_dim_eq)
+
+    # set of IR tuples for converters with amplitude-penalised outputs
+
+    def init_set_IR_dim_neg(m):
+        return ((i, r_i) for (i, r_i) in m.set_IR if r_i in m.set_R_dim_neg[i])
+
+    model.set_IR_dim_neg = pyo.Set(dimen=2, initialize=init_set_IR_dim_neg)
+
+    # set of IR tuples for converters with amplitude-penalised outputs
+
+    def init_set_IR_dim(m):
+        return ((i, r_i) for (i, r_i) in m.set_IR if r_i in m.set_R_dim[i])
+
+    model.set_IR_dim = pyo.Set(dimen=2, initialize=init_set_IR_dim)
+
+    # set of IR tuples for converters with pos. amplitude-constrained outputs
+
+    def init_set_IR_dim_pos(m):
+        return ((i, r_i) for (i, r_i) in m.set_IR if r_i in m.set_R_dim_pos[i])
+
+    model.set_IR_dim_pos = pyo.Set(dimen=2, initialize=init_set_IR_dim_pos)
+
+    # set of IR tuples for converters with externality-inducing outputs
+
+    def init_set_IR_ext(m):
+        return ((i, r_i) for (i, r_i) in m.set_IR if r_i in m.set_R_ext[i])
+
+    model.set_IR_ext = pyo.Set(dimen=2, initialize=init_set_IR_ext)
+
+    # *************************************************************************
+
+    # combined inputs/states/outputs
+
+    # TODO: narrow down these sets if possible
+
+    # set of INN tuples
+
+    def init_set_INN(m):
+        return ((i, n1, n2) for (i, n1) in m.set_IN for n2 in m.set_N[i])
+
+    model.set_INN = pyo.Set(dimen=3, initialize=init_set_INN)
+
+    # set of INM tuples
+
+    def init_set_INM(m):
+        return ((i, n_i, m_i) for (i, n_i) in m.set_IN for m_i in m.set_M[i])
+
+    model.set_INM = pyo.Set(dimen=3, initialize=init_set_INM)
+
+    # set of IRM tuples
+
+    def init_set_IRM(m):
+        return (
+            (i, r_i, m_i) for (i, r_i) in m.set_IR for m_i in m.set_M[i]
+        )  # can be further constrained
+
+    model.set_IRM = pyo.Set(dimen=3, initialize=init_set_IRM)
+    # set of IRN tuples
+
+    def init_set_IRN(m):
+        return (
+            (i, r_i, n_i) for (i, r_i) in m.set_IR for n_i in m.set_N[i]
+        )  # can be further constrained
+
+    model.set_IRN = pyo.Set(dimen=3, initialize=init_set_IRN)
+    
+    # *************************************************************************
+    # *************************************************************************
+
+    # parameters
+
+    # converters
+
+    # externality cost per input unit
+
+    model.param_c_ext_u_imqk = pyo.Param(
+        model.set_IM_ext, model.set_QK, within=pyo.NonNegativeReals, default=0
+    )
+
+    # externality cost per output unit
+
+    model.param_c_ext_y_irqk = pyo.Param(
+        model.set_IR_ext, model.set_QK, within=pyo.NonNegativeReals, default=0
+    )
+
+    # externality cost per state unit
+
+    model.param_c_ext_x_inqk = pyo.Param(
+        model.set_IN_ext, model.set_QK, within=pyo.NonNegativeReals, default=0
+    )
+
+    # unit cost of positive state variations
+
+    model.param_c_pos_var_in = pyo.Param(
+        model.set_IN_pos_var, within=pyo.NonNegativeReals, default=0
+    )
+
+    # unit cost of negative state variations
+
+    model.param_c_neg_var_in = pyo.Param(
+        model.set_IN_neg_var, within=pyo.NonNegativeReals, default=0
+    )
+
+    # unit cost of upper state reference violations
+
+    model.param_c_ref_u_inqk = pyo.Param(
+        model.set_IN_ref_u, model.set_QK, within=pyo.NonNegativeReals, default=0
+    )
+
+    # unit cost of lower state reference violations
+
+    model.param_c_ref_d_inqk = pyo.Param(
+        model.set_IN_ref_d, model.set_QK, within=pyo.NonNegativeReals, default=0
+    )
+
+    # minimum converter cost
+
+    model.param_c_cvt_min_i = pyo.Param(
+        model.set_I_new, within=pyo.NonNegativeReals, default=0
+    )
+
+    # unit (positive) input amplitude cost
+
+    model.param_c_cvt_u_im = pyo.Param(
+        model.set_IM_dim, within=pyo.NonNegativeReals, default=0
+    )
+
+    # unit output amplitude cost
+
+    model.param_c_cvt_y_ir = pyo.Param(
+        model.set_IR_dim, within=pyo.NonNegativeReals, default=0
+    )
+
+    # unit positive state amplitude cost
+
+    model.param_c_cvt_x_pos_in = pyo.Param(
+        model.set_IN_dim_pos, within=pyo.NonNegativeReals, default=0
+    )
+
+    # unit negative state amplitude cost
+
+    model.param_c_cvt_x_neg_in = pyo.Param(
+        model.set_IN_dim_neg, within=pyo.NonNegativeReals, default=0
+    )
+
+    # unit positive output amplitude cost
+
+    model.param_c_cvt_y_pos_ir = pyo.Param(
+        model.set_IR_dim_pos, within=pyo.NonNegativeReals, default=0
+    )
+
+    # unit negative output amplitude cost
+
+    model.param_c_cvt_y_neg_ir = pyo.Param(
+        model.set_IR_dim_neg, within=pyo.NonNegativeReals, default=0
+    )
+
+    # *************************************************************************
+
+    # effect of system inputs on specific network and node pairs
+
+    model.param_a_nw_glimqk = pyo.Param(
+        model.set_GL_not_exp_imp,
+        model.set_IM,
+        model.set_QK,
+        default=0,  # default: no effect
+        within=pyo.Reals,
+    )
+
+    # effect of system outputs on specific network and node pairs
+
+    model.param_a_nw_glirqk = pyo.Param(
+        model.set_GL_not_exp_imp,
+        model.set_IR,
+        model.set_QK,
+        default=0,  # default: no effect
+        within=pyo.Reals,
+    )
+
+    # *************************************************************************
+
+    # inputs
+
+    # upper bounds for (non-binary, non-dimensionable) inputs
+
+    model.param_u_ub_imqk = pyo.Param(
+        model.set_IM_fix, model.set_QK, within=pyo.PositiveReals
+    )
+
+    # maximum input limits
+
+    model.param_u_amp_max_im = pyo.Param(
+        model.set_IM_dim, within=pyo.PositiveReals, default=1
+    )
+
+    # time interval-dependent adjustment coefficients for input limits
+
+    model.param_f_amp_u_imqk = pyo.Param(
+        model.set_IM_dim, model.set_QK, within=pyo.PositiveReals, default=1
+    )
+
+    # *************************************************************************
+
+    # states
+
+    # initial conditions
+
+    model.param_x_inq0 = pyo.Param(model.set_IN, model.set_Q, within=pyo.Reals)
+
+    # fixed upper bounds for state variables
+
+    model.param_x_ub_irqk = pyo.Param(model.set_IN_fix, model.set_QK, within=pyo.Reals)
+
+    # fixed lower bounds for state variables
+
+    model.param_x_lb_irqk = pyo.Param(model.set_IN_fix, model.set_QK, within=pyo.Reals)
+
+    # maximum positive amplitude for states
+
+    model.param_x_amp_pos_max_in = pyo.Param(
+        model.set_IN_dim_pos, within=pyo.PositiveReals
+    )
+
+    # maximum negative amplitude for states
+
+    model.param_x_amp_neg_max_in = pyo.Param(
+        model.set_IN_dim_neg, within=pyo.PositiveReals
+    )
+
+    # adjustment of positive state amplitude limits
+
+    model.param_f_amp_pos_x_inqk = pyo.Param(
+        model.set_IN_dim_pos, model.set_QK, within=pyo.PositiveReals, default=1
+    )
+
+    # adjustment of negative state amplitude limits
+
+    model.param_f_amp_neg_x_inqk = pyo.Param(
+        model.set_IN_dim_neg, model.set_QK, within=pyo.PositiveReals, default=1
+    )
+
+    # state equations: coefficients from C matrix
+
+    model.param_a_eq_x_innqk = pyo.Param(
+        model.set_INN, model.set_QK, default=0, within=pyo.Reals  # default: no effect
+    )
+
+    # state equations: coefficients from D matrix
+
+    model.param_b_eq_x_inmqk = pyo.Param(
+        model.set_INM, model.set_QK, default=0, within=pyo.Reals  # default: no effect
+    )
+
+    # state equations: constant term
+
+    model.param_e_eq_x_inqk = pyo.Param(
+        model.set_IN, model.set_QK, default=0, within=pyo.Reals  # default: no effect
+    )
+
+    # *************************************************************************
+
+    # outputs
+
+    # fixed upper bounds for output variables
+
+    model.param_y_ub_irqk = pyo.Param(model.set_IR_fix, model.set_QK, within=pyo.Reals)
+
+    # fixed lower bounds for output variables
+
+    model.param_y_lb_irqk = pyo.Param(model.set_IR_fix, model.set_QK, within=pyo.Reals)
+
+    # adjustment of positive output amplitude limits
+
+    model.param_f_amp_y_pos_irqk = pyo.Param(
+        model.set_IR_dim_pos, model.set_QK, within=pyo.PositiveReals, default=1
+    )
+
+    # adjustment of negative output amplitude limits
+
+    model.param_f_amp_y_neg_irqk = pyo.Param(
+        model.set_IR_dim_neg, model.set_QK, within=pyo.PositiveReals, default=1
+    )
+
+    # maximum positive amplitude limit for outputs
+
+    model.param_y_amp_pos_max_ir = pyo.Param(
+        model.set_IR_dim_pos, within=pyo.PositiveReals
+    )
+
+    # maximum negative amplitude limit for outputs
+
+    model.param_y_amp_neg_max_ir = pyo.Param(
+        model.set_IR_dim_neg, within=pyo.PositiveReals
+    )
+
+    # output equation coefficients from C matrix
+
+    model.param_c_eq_y_irnqk = pyo.Param(
+        model.set_IRN, model.set_QK, default=0, within=pyo.Reals  # default: no effect
+    )
+
+    # output equation coefficients from D matrix
+
+    model.param_d_eq_y_irmqk = pyo.Param(
+        model.set_IRM, model.set_QK, default=0, within=pyo.Reals  # default: no effect
+    )
+
+    # output equation constant
+
+    model.param_e_eq_y_irqk = pyo.Param(
+        model.set_IR, model.set_QK, default=0, within=pyo.Reals  # default: no effect
+    )
+
+    # *************************************************************************
+    # *************************************************************************
+    # *************************************************************************
+    # *************************************************************************
+
+    # variables
+
+    # *************************************************************************
+    # *************************************************************************
+
+    # capex for installing individual converters
+
+    model.var_capex_cvt_i = pyo.Var(model.set_I_new, within=pyo.NonNegativeReals)
+
+    # *************************************************************************
+
+    # converters
+
+    # decision to install converter i
+
+    model.var_cvt_inv_i = pyo.Var(model.set_I_new, within=pyo.Binary)
+
+    # inputs
+
+    # input variables
+
+    def bounds_var_u_imqk(m, i, m_i, q, k):
+        if (i, m_i) in m.param_u_ub_imqk:
+            # predefined limit
+            return (0, m.param_u_ub_imqk[(i, m_i, q, k)])
+        else:
+            # dynamic limit (set elsewhere)
+            return (0, None)
+
+    def domain_var_u_imqk(m, i, m_i, q, k):
+        try:
+            if m_i in m.set_M_bin[i]:
+                return pyo.Binary  # binary: {0,1}
+            else:
+                return pyo.NonNegativeReals  # nonnegative real: [0,inf]
+        except KeyError:
+            return pyo.NonNegativeReals  # nonnegative real: [0,inf]
+
+    model.var_u_imqk = pyo.Var(
+        model.set_IM,
+        model.set_QK,
+        domain=domain_var_u_imqk,
+        # within=pyo.NonNegativeReals,
+        bounds=bounds_var_u_imqk,
+    )
+
+    # input amplitude variables (only one per sign is needed, as vars. are nnr)
+
+    model.var_u_amp_im = pyo.Var(model.set_IM_dim, within=pyo.NonNegativeReals)
+
+    # *************************************************************************
+
+    # outputs
+
+    # output variables
+
+    def bounds_var_y_irqk(m, i, r, q, k):
+        if r in m.set_R_fix:
+            # predefined limit
+            return (m.param_u_lb_irqk[(i, r, q, k)], m.param_u_ub_irqk[(i, r, q, k)])
+        else:
+            # do not enforce any limits
+            return (None, None)
+
+    # def domain_var_y_irqk(m, i, r, k):
+    #     try:
+    #         if m_i in m.set_M_bin[i]:
+    #             return pyo.Binary # binary: {0,1}
+    #         else:
+    #             return pyo.NonNegativeReals # nonnegative real: [0,inf]
+    #     except KeyError:
+    #         return pyo.NonNegativeReals # nonnegative real: [0,inf]
+
+    model.var_y_irqk = pyo.Var(
+        model.set_IR, model.set_QK, bounds=bounds_var_y_irqk, within=pyo.Reals
+    )
+
+    # positive output amplitudes
+
+    model.var_y_amp_pos_ir = pyo.Var(model.set_IR_dim_pos, within=pyo.Reals)
+
+    # output amplitudes
+
+    model.var_y_amp_neg_ir = pyo.Var(model.set_IR_dim_neg, within=pyo.Reals)
+
+    # *************************************************************************
+
+    # states
+
+    # state variables
+
+    model.var_x_inqk = pyo.Var(model.set_IN, model.set_QK, within=pyo.Reals)
+
+    # positive amplitude variables
+
+    model.var_x_amp_pos_in = pyo.Var(model.set_IN_dim_pos, within=pyo.NonNegativeReals)
+
+    # negative amplitude variables
+
+    model.var_x_amp_neg_in = pyo.Var(model.set_IN_dim_neg, within=pyo.NonNegativeReals)
+
+    # positive state variation
+
+    model.var_delta_x_pos_var_in = pyo.Var(
+        model.set_IN_pos_var, within=pyo.NonNegativeReals
+    )
+
+    # negative state variation
+
+    model.var_delta_x_neg_var_in = pyo.Var(
+        model.set_IN_neg_var, within=pyo.NonNegativeReals
+    )
+
+    # positive reference state violation
+
+    model.var_delta_x_ref_u_inqk = pyo.Var(
+        model.set_IN_ref_u, model.set_QK, within=pyo.NonNegativeReals
+    )
+
+    # negative reference state violation
+
+    model.var_delta_x_ref_d_inqk = pyo.Var(
+        model.set_IN_ref_d, model.set_QK, within=pyo.NonNegativeReals
+    )
+
+    # *************************************************************************
+    # *************************************************************************
+
+    # objective function
+
+
+    # capex for converters
+
+    def rule_capex_converter(m, i):
+        return (
+            m.var_cvt_inv_i[i] * m.param_c_cvt_min_i[i]
+            + sum(
+                m.var_u_amp_im[(i, m_i)] * m.param_c_cvt_u_im[(i, m_i)]
+                for m_i in m.set_M_dim_i[i]
+            )
+            + sum(
+                m.var_x_amp_pos_in[(i, n)] * m.param_c_cvt_x_pos_in[(i, n)]
+                for n in m.set_N_dim_pos[i]
+            )
+            + sum(
+                m.var_x_amp_neg_in[(i, n)] * m.param_c_cvt_x_neg_in[(i, n)]
+                for n in m.set_N_dim_neg[i]
+            )
+            + sum(
+                m.var_y_amp_pos_ir[(i, r)] * m.param_c_cvt_y_pos_ir[(i, r)]
+                for r in m.set_N_dim_pos[i]
+            )
+            + sum(
+                m.var_y_amp_neg_ir[(i, r)] * m.param_c_cvt_y_neg_ir[(i, r)]
+                for r in m.set_N_dim_neg[i]
+            )
+            <= m.var_capex_cvt_i[i]
+        )
+
+    model.constr_capex_system = pyo.Constraint(
+        model.set_I_new, rule=rule_capex_converter
+    )
+
+    # *************************************************************************
+    # *************************************************************************
+
+    # converters
+
+    # *************************************************************************
+    # *************************************************************************
+
+    # input signal limits for dimensionable inputs
+
+    def rule_constr_u_limit_dim(m, i, m_i, q, k):
+        return (
+            m.var_u_imqk[(i, m_i, q, k)]
+            <= m.var_u_amp_im[(i, m_i)] * m.param_f_amp_u_imqk[(i, m_i, q, k)]
+        )
+
+    model.constr_u_limit_dim = pyo.Constraint(
+        model.set_IM_dim, model.set_QK, rule=rule_constr_u_limit_dim
+    )
+
+    # nominal input amplitude limit for dimensionable inputs
+
+    def rule_constr_u_amp_ub(m, i, m_i):
+        return (
+            m.var_u_amp_im[(i, m_i)]
+            <= m.var_cvt_inv_i[i] * m.param_u_amp_max_im[(i, m_i)]
+        )
+
+    model.constr_u_amp_ub = pyo.Constraint(model.set_IM_dim, rule=rule_constr_u_amp_ub)
+
+    # fixed upper limits
+
+    def rule_constr_u_fix_limits(m, i, m_i, q, k):
+        # if we need to know the lim input signal (e.g., for the obj. func.)
+
+        if i in m.set_I_new:
+            # new converter
+
+            return (
+                m.var_u_imqk[(i, m_i, q, k)]
+                <= m.param_u_ub_imqk[(i, m_i, q, k)] * m.var_cvt_inv_i[i]
+            )
+
+            return m.var_u_imqk[(i, m_i, q, k)] <= m.var_cvt_inv_i[i]
+
+        else:
+            # pre-existing
+
+            return m.var_u_imqk[(i, m_i, q, k)] <= m.param_u_ub_imqk[(i, m_i, q, k)]
+
+    model.constr_u_fix_limits = pyo.Constraint(
+        model.set_IM_fix, model.set_QK, rule=rule_constr_u_fix_limits
+    )
+
+    # input limits for binary inputs
+
+    def rule_constr_u_bin_limits(m, i, m_i, q, k):
+        if i in m.set_I_new:
+            # binary variables
+
+            return m.var_u_imqk[(i, m_i, q, k)] <= m.var_cvt_inv_i[i]
+
+        else:
+            return pyo.Constraint.Skip
+
+    model.constr_u_bin_limits = pyo.Constraint(
+        model.set_IM_bin, model.set_QK, rule=rule_constr_u_bin_limits
+    )
+
+    # *************************************************************************
+
+    # outputs
+
+    # output equations
+
+    def rule_constr_output_equations(m, i, r, q, k):
+        return (
+            m.var_y_irqk[(i, r, k)]
+            == sum(
+                m.param_c_eq_y_irnqk[(i, r, n_i, q, k)] * m.var_x_inqk[(i, n_i, q, k)]
+                for n_i in m.set_N[i]
+            )
+            + sum(
+                m.param_d_eq_y_irmqk[(i, r, m_i, q, k)] * m.var_u_imqk[(i, m_i, q, k)]
+                for m_i in m.set_M[i]
+            )
+            + m.param_e_eq_y_irqk[(i, r, q, k)]
+        )
+
+    model.constr_output_equations = pyo.Constraint(
+        model.set_IR, model.set_QK, rule=rule_constr_output_equations
+    )
+
+    # positive amplitude limit for output variables
+
+    def rule_constr_y_vars_have_pos_amp_limits(m, i, r, q, k):
+        return m.var_y_irqk[(i, r, q, k)] <= (
+            m.var_y_amp_pos_ir[(i, r)] * m.param_f_amp_y_pos_irqk[(i, r, q, k)]
+        )
+
+    model.constr_y_vars_have_pos_amp_limits = pyo.Constraint(
+        model.set_IR_dim_pos, model.set_QK, rule=rule_constr_y_vars_have_pos_amp_limits
+    )
+
+    # negative amplitude limit for output variables
+
+    def rule_constr_y_vars_have_neg_amp_limits(m, i, r, q, k):
+        return m.var_y_irqk[(i, r, q, k)] >= (
+            -m.var_y_amp_neg_ir[(i, r)] * m.param_f_amp_y_neg_irqk[(i, r, q, k)]
+        )
+
+    model.constr_y_vars_have_neg_amp_limits = pyo.Constraint(
+        model.set_IR_dim_neg, model.set_QK, rule=rule_constr_y_vars_have_neg_amp_limits
+    )
+
+    # positive amplitude limit must be zero unless the system is installed
+
+    def rule_constr_y_amp_pos_zero_if_cvt_not_selected(m, i, r):
+        return m.var_y_amp_pos_ir[(i, r)] <= (
+            m.var_cvt_inv_i[i] * m.param_y_amp_pos_ir[(i, r)]
+        )
+
+    model.constr_y_amp_pos_zero_if_cvt_not_newected = pyo.Constraint(
+        model.set_IR_dim_pos, rule=rule_constr_y_amp_pos_zero_if_cvt_not_selected
+    )
+
+    # negative amplitude limit must be zero unless the system is installed
+
+    def rule_constr_y_amp_neg_zero_if_cvt_not_selected(m, i, r):
+        return m.var_y_amp_neg_ir[(i, r)] <= (
+            m.var_cvt_inv_i[i] * m.param_y_amp_neg_ir[(i, r)]
+        )
+
+    model.constr_y_amp_neg_zero_if_cvt_not_selected = pyo.Constraint(
+        model.set_IR_dim_neg, rule=rule_constr_y_amp_neg_zero_if_cvt_not_selected
+    )
+
+    # the positive and negative amplitudes must match
+
+    def rule_constr_y_amp_pos_neg_match(m, i, r):
+        return m.var_y_amp_pos_ir[(i, r)] == m.var_y_amp_neg_ir[(i, r)]
+
+    model.constr_y_amp_pos_neg_match = pyo.Constraint(
+        model.set_IR_dim_eq, rule=rule_constr_y_amp_pos_neg_match
+    )
+
+    # *************************************************************************
+
+    # states
+
+    def rule_constr_state_equations(m, i, n, q, k):
+        return (
+            m.var_x_inqk[(i, n, q, k)]
+            == sum(
+                m.param_a_eq_x_innqk[(i, n, n_star, q, k)]
+                * (
+                    m.var_x_inqk[(i, n_star, q, k - 1)]
+                    if k != 0
+                    else m.param_x_inq0[(i, n, q)]
+                )
+                for n_star in m.set_N[i]
+            )
+            + sum(
+                m.param_b_eq_x_inmqk[(i, n, m_i, q, k)] * m.var_u_imqk[(i, m_i, q, k)]
+                for m_i in m.set_M[i]
+            )
+            + m.param_e_eq_x_inqk[(i, n, q, k)]
+        )
+
+    model.constr_state_equations = pyo.Constraint(
+        model.set_IN, model.set_QK, rule=rule_constr_state_equations
+    )
+
+    # positive amplitude limit for state variables
+
+    def rule_constr_x_vars_have_pos_amp_limits(m, i, n, q, k):
+        return m.var_x_inqk[(i, n, q, k)] <= (
+            m.var_x_amp_pos_in[(i, n)] * m.param_f_amp_x_pos_inqk[(i, n, q, k)]
+        )
+
+    model.constr_x_vars_have_pos_amp_limits = pyo.Constraint(
+        model.set_IN_dim_pos, model.set_QK, rule=rule_constr_x_vars_have_pos_amp_limits
+    )
+
+    # negative amplitude limit for state variables
+
+    def rule_constr_x_vars_have_neg_amp_limits(m, i, n, q, k):
+        return m.var_x_inqk[(i, n, q, k)] >= (
+            -m.var_y_amp_neg_in[(i, n)] * m.param_f_amp_x_neg_inqk[(i, n, q, k)]
+        )
+
+    model.constr_x_vars_have_neg_amp_limits = pyo.Constraint(
+        model.set_IN_dim_neg, model.set_QK, rule=rule_constr_x_vars_have_neg_amp_limits
+    )
+
+    # positive amplitude limit must be zero unless the system is installed
+
+    def rule_constr_x_amp_pos_zero_if_cvt_not_selected(m, i, n):
+        return m.var_x_amp_pos_in[(i, n)] <= (
+            m.var_cvt_inv_i[i] * m.param_x_amp_pos_in[(i, n)]
+        )
+
+    model.constr_x_amp_pos_zero_if_cvt_not_selected = pyo.Constraint(
+        model.set_IN_dim_pos, rule=rule_constr_x_amp_pos_zero_if_cvt_not_selected
+    )
+
+    # negative amplitude limit must be zero unless the system is installed
+
+    def rule_constr_x_amp_neg_zero_if_cvt_not_selected(m, i, n):
+        return m.var_x_amp_neg_in[(i, n)] <= (
+            m.var_cvt_inv_i[i] * m.param_x_amp_neg_in[(i, n)]
+        )
+
+    model.constr_x_amp_neg_zero_if_cvt_not_selected = pyo.Constraint(
+        model.set_IN_dim_neg, rule=rule_constr_x_amp_neg_zero_if_cvt_not_selected
+    )
+
+    # the positive and negative amplitudes must match
+
+    def rule_constr_x_amp_pos_neg_match(m, i, n):
+        return m.var_x_amp_pos_in[(i, n)] == m.var_x_amp_neg_in[(i, n)]
+
+    model.constr_x_amp_pos_neg_match = pyo.Constraint(
+        model.set_IN_dim_eq, rule=rule_constr_x_amp_pos_neg_match
+    )
+
+    # *************************************************************************
+    # *************************************************************************
+
+    return model
+
+    # *************************************************************************
+    # *************************************************************************
+
+
+# *****************************************************************************
+# *****************************************************************************
+# *****************************************************************************
+# *****************************************************************************
+# *****************************************************************************
diff --git a/src/topupopt/problems/esipp/blocks/networks.py b/src/topupopt/problems/esipp/blocks/networks.py
new file mode 100644
index 0000000..1a3823b
--- /dev/null
+++ b/src/topupopt/problems/esipp/blocks/networks.py
@@ -0,0 +1,746 @@
+# imports
+
+import pyomo.environ as pyo
+
+from math import isfinite, inf
+
+# *****************************************************************************
+# *****************************************************************************
+
+
+def add_network_restrictions(
+        model: pyo.AbstractModel,
+        enable_default_values: bool = True,
+        enable_validation: bool = True,
+        enable_initialisation: bool = True,
+):
+    
+    # *************************************************************************
+    # *************************************************************************
+    
+    model.set_L_max_in_g = pyo.Set(
+        model.set_G, within=model.set_L
+    )  # should inherently exclude import nodes
+
+    model.set_L_max_out_g = pyo.Set(
+        model.set_G, within=model.set_L
+    )  # should inherently exclude export nodes
+
+    # maximum number of arcs per node pair
+
+    model.param_max_number_parallel_arcs = pyo.Param(
+        model.set_GLL,
+        # within=pyo.PositiveIntegers,
+        within=pyo.PositiveReals,
+        default=inf,
+    )
+
+    def init_set_GLL_arc_max(m):
+        return (
+            (g, l1, l2)
+            for (g, l1, l2) in m.param_max_number_parallel_arcs
+            if isfinite(m.param_max_number_parallel_arcs[(g, l1, l2)])
+        )
+
+    model.set_GLL_arc_max = pyo.Set(
+        dimen=3, within=model.set_GLL, initialize=init_set_GLL_arc_max
+    )
+
+    # *************************************************************************
+    # *************************************************************************
+
+    # limit number of directed arcs per direction
+
+    def rule_constr_limited_parallel_arcs_per_direction(m, g, l1, l2):
+        # cases:
+        # 1) the number of options is lower than or equal to the limit (skip)
+        # 2) the number of preexisting and new mandatory arcs exceeds
+        # the limit (infeasible: pyo.Constraint.Infeasible)
+        # 3) all other cases (constraint)
+
+        # number of preexisting arcs going from l1 to l2
+
+        number_arcs_pre_nom = (
+            len(m.set_J_pre[(g, l1, l2)]) if (g, l1, l2) in m.set_J_pre else 0
+        )
+
+        number_arcs_pre_rev = (
+            sum(1 for j in m.set_J_pre[(g, l2, l1)] if j in m.set_J_und[(g, l2, l1)])
+            if (g, l2, l1) in m.set_J_pre
+            else 0
+        )
+
+        # number of mandatory arcs going from l1 to l2
+
+        number_arcs_mdt_nom = (
+            len(m.set_J_mdt[(g, l1, l2)]) if (g, l1, l2) in m.set_J_mdt else 0
+        )
+
+        number_arcs_mdt_rev = (
+            sum(1 for j in m.set_J_mdt[(g, l2, l1)] if j in m.set_J_und[(g, l2, l1)])
+            if (g, l2, l1) in m.set_J_mdt
+            else 0
+        )
+
+        # number of optional arcs going from l1 to l2
+
+        number_arcs_opt_nom = (
+            sum(
+                1
+                for j in m.set_J[(g, l1, l2)]
+                if j not in m.set_J_pre[(g, l1, l2)]
+                if j not in m.set_J_mdt[(g, l1, l2)]
+            )
+            if (g, l1, l2) in m.set_J
+            else 0
+        )
+
+        number_arcs_opt_rev = (
+            sum(
+                1
+                for j in m.set_J[(g, l2, l1)]
+                if j not in m.set_J_pre[(g, l2, l1)]
+                if j not in m.set_J_mdt[(g, l2, l1)]
+                if j in m.set_J_und[(g, l2, l1)]
+            )
+            if (g, l2, l1) in m.set_J
+            else 0
+        )
+
+        # build the constraints
+
+        if (
+            number_arcs_mdt_nom
+            + number_arcs_mdt_rev
+            + number_arcs_pre_nom
+            + number_arcs_pre_rev
+            > m.param_max_number_parallel_arcs[(g, l1, l2)]
+        ):
+            # the number of unavoidable arcs already exceeds the limit
+            return pyo.Constraint.Infeasible
+
+        elif (
+            number_arcs_opt_nom
+            + number_arcs_opt_rev
+            + number_arcs_mdt_nom
+            + number_arcs_mdt_rev
+            + number_arcs_pre_nom
+            + number_arcs_pre_rev
+            > m.param_max_number_parallel_arcs[(g, l1, l2)]
+        ):
+            # the number of potential arcs exceeds the limit: cannot be skipped
+            return (
+                # preexisting arcs
+                number_arcs_pre_nom + number_arcs_pre_rev +
+                # mandatory arcs
+                number_arcs_mdt_nom + number_arcs_mdt_rev +
+                # arcs within an (optional) group that uses interfaces
+                sum(
+                    (
+                        sum(
+                            1
+                            for j in m.set_J_col[(g, l1, l2)]
+                            if (g, l1, l2, j) in m.set_GLLJ_col_t[t]
+                        )
+                        if (g, l1, l2) in m.set_J_col
+                        else 0
+                        + sum(
+                            1
+                            for j in m.set_J_col[(g, l2, l1)]
+                            if j in m.set_J_und[(g, l2, l1)]
+                            if (g, l2, l1, j) in m.set_GLLJ_col_t[t]
+                        )
+                        if ((g, l2, l1) in m.set_J_col and (g, l2, l1) in m.set_J_und)
+                        else 0
+                    )
+                    * m.var_xi_arc_inv_t[t]
+                    for t in m.set_T_int
+                )
+                +
+                # arcs within an (optional) group that does not use interfaces
+                sum(
+                    (
+                        sum(
+                            1
+                            for j in m.set_J_col[(g, l1, l2)]
+                            if (g, l1, l2, j) in m.set_GLLJ_col_t[t]
+                        )
+                        if (g, l1, l2) in m.set_J_col
+                        else 0
+                        + sum(
+                            1
+                            for j in m.set_J_col[(g, l2, l1)]
+                            if j in m.set_J_und[(g, l2, l1)]
+                            if (g, l2, l1, j) in m.set_GLLJ_col_t[t]
+                        )
+                        if ((g, l2, l1) in m.set_J_col and (g, l2, l1) in m.set_J_und)
+                        else 0
+                    )
+                    * sum(m.var_delta_arc_inv_th[(t, h)] for h in m.set_H_t[t])
+                    for t in m.set_T  # new
+                    if t not in m.set_T_mdt  # optional
+                    if t not in m.set_T_int  # not interfaced
+                )
+                +
+                # optional individual arcs using interfaces, nominal direction
+                sum(
+                    m.var_xi_arc_inv_gllj[(g, l1, l2, j)]
+                    for j in m.set_J_int[(g, l1, l2)]  # interfaced
+                    if j not in m.set_J_col[(g, l1, l2)]  # individual
+                )
+                if (g, l1, l2) in m.set_J_int
+                else 0 +
+                # optional individual arcs using interfaces, reverse direction
+                sum(
+                    m.var_xi_arc_inv_gllj[(g, l2, l1, j)]
+                    for j in m.set_J_int[(g, l2, l1)]  # interfaced
+                    if j in m.set_J_und[(g, l2, l1)]  # undirected
+                    if j not in m.set_J_col[(g, l1, l2)]  # individual
+                )
+                if ((g, l2, l1) in m.set_J_int and (g, l2, l1) in m.set_J_und)
+                else 0 +
+                # optional individual arcs not using interfaces, nominal dir.
+                sum(
+                    sum(
+                        m.var_delta_arc_inv_glljh[(g, l1, l2, j, h)]
+                        for h in m.set_H_gllj[(g, l1, l2, j)]
+                    )
+                    for j in m.set_J[(g, l1, l2)]
+                    if j not in m.set_J_pre[(g, l1, l2)]  # not preexisting
+                    if j not in m.set_J_mdt[(g, l1, l2)]  # not mandatory
+                    if j not in m.set_J_int[(g, l1, l2)]  # not interfaced
+                    if j not in m.set_J_col[(g, l1, l2)]  # individual
+                )
+                if (g, l1, l2) in m.set_J
+                else 0 +
+                # optional individual arcs not using interfaces, reverse dir.
+                sum(
+                    sum(
+                        m.var_delta_arc_inv_glljh[(g, l2, l1, j, h)]
+                        for h in m.set_H_gllj[(g, l2, l1, j)]
+                    )
+                    for j in m.set_J_opt[(g, l2, l1)]
+                    if j in m.set_J_und[(g, l2, l1)]
+                    if j not in m.set_J_pre[(g, l2, l1)]  # not preexisting
+                    if j not in m.set_J_mdt[(g, l2, l1)]  # not mandatory
+                    if j not in m.set_J_int[(g, l2, l1)]  # not interfaced
+                    if j not in m.set_J_col[(g, l2, l1)]  # individual
+                )
+                if (g, l2, l1) in m.set_J
+                else 0 <= m.param_max_number_parallel_arcs[(g, l1, l2)]
+            )
+
+        else:  # the number of options is lower than or equal to the limit: skip
+            return pyo.Constraint.Skip
+
+    model.constr_limited_parallel_arcs_per_direction = pyo.Constraint(
+        model.set_GLL_arc_max, rule=rule_constr_limited_parallel_arcs_per_direction
+    )
+
+    # *************************************************************************
+    # *************************************************************************
+
+    # there can only one incoming arc at most, if there are no outgoing arcs
+
+    def rule_constr_max_incoming_directed_arcs(m, g, l):
+        # check if the node is not among those subject to a limited number of incoming arcs
+        if l not in m.set_L_max_in_g[g]:
+            # it is not, skip this constraint
+            return pyo.Constraint.Skip
+
+        # max number of directed incoming arcs
+        n_max_dir_in = sum(
+            sum(
+                1
+                for j in m.set_J[(g, l_line, l)]
+                if j not in m.set_J_und[(g, l_line, l)]
+            )  # directed
+            for l_line in m.set_L[g] # for every node
+            if l_line != l # cannot be the same node
+            # if l_line not in m.set_L_imp[g] # why?
+            if (g, l_line, l) in m.set_J
+        )
+        
+        # check the maximum number of incoming arcs
+        if n_max_dir_in <= 1:
+            # there can only be one incoming arc at most: redundant constraint
+            return pyo.Constraint.Skip
+        else:  # more than one incoming arc is possible
+        
+            # number of (new) incoming directed arcs in a group
+            b_max_in_gl = 0
+
+            # the big m
+
+            M_gl = n_max_dir_in - 1  # has to be positive since n_max_dir_in > 1
+            # TODO: put parenthesis to avoid funny results
+            temp_constr = (
+                sum(
+                    # *********************************************************
+                    # interfaced groups
+                    sum(
+                        sum(
+                            1
+                            for j in m.set_J_col[(g, l_circ, l)]  # part of group
+                            if j not in m.set_J_und[(g, l_circ, l)]  # directed
+                            if (g, l_circ, l, j) in m.set_GLLJ_col_t[t]
+                        )
+                        * m.var_xi_arc_inv_t[t]  # in t
+                        for t in m.set_T_int
+                    )
+                    +
+                    # *********************************************************
+                    # optional non-interfaced groups
+                    sum(
+                        sum(
+                            sum(
+                                1
+                                for j in m.set_J_col[(g, l_circ, l)]  # part of group
+                                if j not in m.set_J_und[(g, l_circ, l)]  # directed
+                                if (g, l_circ, l, j) in m.set_GLLJ_col_t[t]
+                            )
+                            * m.var_delta_arc_inv_th[(t, h)]
+                            for h in m.set_H_t[t]
+                        )
+                        for t in m.set_T
+                        if t not in m.set_T_mdt  # optional
+                        if t not in m.set_T_int  # not interfaced
+                    )
+                    +
+                    # *********************************************************
+                    # interfaced arcs
+                    (sum(
+                        m.var_xi_arc_inv_gllj[(g, l_circ, l, j_circ)]
+                        for j_circ in m.set_J[(g, l_circ, l)]
+                        if j_circ not in m.set_J_und[(g, l_circ, l)]  # directed
+                        if j_circ in m.set_J_int[(g, l_circ, l)]  # interfaced
+                        if j_circ not in m.set_J_col[(g, l_circ, l)]  # individual
+                    )
+                    if (g, l_circ, l) in m.set_J
+                    else 0) +
+                    # *********************************************************
+                    # optional non-interfaced arcs
+                    (sum(
+                        sum(
+                            m.var_delta_arc_inv_glljh[(g, l_circ, l, j_dot, h_dot)]
+                            for h_dot in m.set_H_gllj[(g, l_circ, l, j_dot)]
+                        )
+                        for j_dot in m.set_J[(g, l_circ, l)]
+                        if j_dot not in m.set_J_und[(g, l_circ, l)]  # directed
+                        if j_dot not in m.set_J_int[(g, l_circ, l)]  # not interfaced
+                        if j_dot not in m.set_J_col[(g, l_circ, l)]  # individual
+                        if j_dot not in m.set_J_pre[(g, l_circ, l)]  # new
+                        if j_dot not in m.set_J_mdt[(g, l_circ, l)]  # optional
+                    )
+                    if (g, l_circ, l) in m.set_J
+                    else 0) +
+                    # *********************************************************
+                    # preexisting directed arcs
+                    (sum(
+                        1
+                        for j_pre_dir in m.set_J_pre[(g, l_circ, l)]  # preexisting
+                        if j_pre_dir not in m.set_J_und[(g, l_circ, l)]  # directed
+                    )
+                    if (g, l_circ, l) in m.set_J_pre
+                    else 0) +
+                    # *********************************************************
+                    # mandatory directed arcs
+                    (sum(
+                        1
+                        for j_mdt_dir in m.set_J_mdt[(g, l_circ, l)]
+                        if j_mdt_dir not in m.set_J_und[(g, l_circ, l)]  # directed
+                    )
+                    if (g, l_circ, l) in m.set_J_mdt
+                    else 0)
+                    # *********************************************************
+                    for l_circ in m.set_L[g]
+                    if l_circ not in m.set_L_exp[g]
+                    if l_circ != l
+                )
+                <= 1  # +
+                # M_gl*sum(
+                #     # *********************************************************
+                #     # outgoing arcs in interfaced groups, nominal direction
+                #     sum(sum(1
+                #             for j in m.set_J_col[(g,l,l_diamond)]
+                #             #if j in m.set_J_int[(g,l,l_diamond)]
+                #             if (g,l,l_diamond,j) in m.set_GLLJ_col_t[t]
+                #             )*m.var_xi_arc_inv_t[t]
+                #         for t in m.set_T_int
+                #         ) if (g,l,l_diamond) in m.set_J_col else 0
+                #     +
+                #     # outgoing arcs in interfaced groups, reverse direction
+                #     sum(sum(1
+                #             for j in m.set_J_col[(g,l_diamond,l)]
+                #             #if j in m.set_J_int[(g,l_diamond,l)]
+                #             if j in m.set_J_und[(g,l_diamond,l)]
+                #             if (g,l_diamond,l,j) in m.set_GLLJ_col_t[t]
+                #             )*m.var_xi_arc_inv_t[t]
+                #         for t in m.set_T_int
+                #         ) if (g,l_diamond,l) in m.set_J_col else 0
+                #     +
+                #     # *********************************************************
+                #     # TODO: outgoing arcs in non-interfaced optional groups, nominal
+                #     sum(sum(1
+                #             for j in m.set_J_col[(g,l,l_diamond)]
+                #             #if j in m.set_J_int[(g,l,l_diamond)]
+                #             if (g,l,l_diamond,j) in m.set_GLLJ_col_t[t]
+                #             )*sum(
+                #                 m.var_delta_arc_inv_th[(t,h)]
+                #                 for h in m.set_H_t[t]
+                #                 )
+                #         for t in m.set_T
+                #         if t not in m.set_T_mdt
+                #         if t not in m.set_T_int
+                #         ) if (g,l,l_diamond) in m.set_J_col else 0
+                #     +
+                #     # TODO: outgoing arcs in non-interfaced optional groups, reverse
+                #     sum(sum(1
+                #             for j in m.set_J_col[(g,l_diamond,l)]
+                #             #if j in m.set_J_int[(g,l_diamond,l)]
+                #             if j in m.set_J_und[(g,l_diamond,l)]
+                #             if (g,l_diamond,l,j) in m.set_GLLJ_col_t[t]
+                #             )*sum(
+                #                 m.var_delta_arc_inv_th[(t,h)]
+                #                 for h in m.set_H_t[t]
+                #                 )
+                #         for t in m.set_T
+                #         if t not in m.set_T_mdt
+                #         if t not in m.set_T_int
+                #         ) if (g,l_diamond,l) in m.set_J_col else 0
+                #     +
+                #     # *********************************************************
+                #     # interfaced individual outgoing arcs, nominal direction
+                #     sum(m.var_xi_arc_inv_gllj[(g,l,l_diamond,j)]
+                #         for j in m.set_J_int[(g,l,l_diamond)] # interfaced
+                #         if j not in m.set_J_col[(g,l,l_diamond)] # individual
+                #         ) if (g,l,l_diamond) in m.set_J_int else 0
+                #     +
+                #     # *********************************************************
+                #     # interfaced individual undirected arcs, reverse direction
+                #     sum(m.var_xi_arc_inv_gllj[(g,l,l_diamond,j)]
+                #         for j in m.set_J_und[(g,l_diamond,l)] # undirected
+                #         if j in m.set_J_int[(g,l_diamond,l)] # interfaced
+                #         if j not in m.set_J_col[(g,l_diamond,l)] # individual
+                #         ) if (g,l_diamond,l) in m.set_J_und else 0
+                #     +
+                #     # *********************************************************
+                #     # outgoing non-interfaced individual optional arcs
+                #     sum(
+                #         sum(m.var_delta_arc_inv_glljh[(g,l,l_diamond,j,h)]
+                #             for h in m.set_H_gllj[(g,l,l_diamond,j)])
+                #         for j in m.set_J[(g,l,l_diamond)]
+                #         if j not in m.set_J_col[(g,l,l_diamond)] # individual
+                #         if j not in m.set_J_mdt[(g,l,l_diamond)] # optional
+                #         if j not in m.set_J_int[(g,l,l_diamond)] # interfaced
+                #         ) if (g,l,l_diamond) in m.set_J else 0
+                #     +
+                #     # *********************************************************
+                #     # individual non-interfaced undirected arcs, reverse dir.
+                #     sum(
+                #         sum(m.var_delta_arc_inv_glljh[(g,l_diamond,l,j,h)]
+                #             for h in m.set_H_gllj[(g,l_diamond,l,j)])
+                #         for j in m.set_J_und[(g,l_diamond,l)] # undirected
+                #         if j not in m.set_J_col[(g,l_diamond,l)] # individual
+                #         if j not in m.set_J_mdt[(g,l_diamond,l)] # optional
+                #         if j not in m.set_J_int[(g,l_diamond,l)] # interfaced
+                #         ) if (g,l_diamond,l) in m.set_J_und else 0
+                #     +
+                #     # *********************************************************
+                #     # preselected outgonig arcs, nominal direction
+                #     len(m.set_J_pre[(g,l,l_diamond)]
+                #         ) if (g,l,l_diamond) in m.set_J_pre else 0
+                #     +
+                #     # *********************************************************
+                #     # mandatory outgoing arcs, nominal direction
+                #     len(m.set_J_mdt[(g,l,l_diamond)]
+                #         ) if (g,l,l_diamond) in m.set_J_mdt else 0
+                #     +
+                #     # *********************************************************
+                #     # undirected preselected arcs, reverse direction
+                #     sum(1
+                #         for j in m.set_J_pre[(g,l_diamond,l)]
+                #         if j in m.set_J_und[(g,l_diamond,l)]
+                #         ) if (g,l_diamond,l) in m.set_J_pre else 0
+                #     +
+                #     # *********************************************************
+                #     # undirected mandatory arcs, reverse direction
+                #     sum(1
+                #         for j in m.set_J_mdt[(g,l_diamond,l)]
+                #         if j in m.set_J_und[(g,l_diamond,l)]
+                #         ) if (g,l_diamond,l) in m.set_J_mdt else 0
+                #     # *********************************************************
+                #     for l_diamond in m.set_L[g]
+                #     if l_diamond not in m.set_L_imp[g]
+                #     if l_diamond != l
+                #     )
+            )
+            if type(temp_constr) == bool:
+                # trivial outcome
+                return pyo.Constraint.Feasible if temp_constr else pyo.Constraint.Infeasible
+            else:
+                # constraint is relevant
+                return temp_constr
+
+    model.constr_max_incoming_directed_arcs = pyo.Constraint(
+        model.set_GL, rule=rule_constr_max_incoming_directed_arcs
+    )
+
+    # *************************************************************************
+    # *************************************************************************
+
+    def rule_constr_max_outgoing_directed_arcs(m, g, l):
+        # check if the node is not among those subject to a limited number of outgoing arcs
+        if l not in m.set_L_max_out_g[g]:
+            # it is not, skip this constraint
+            return pyo.Constraint.Skip
+
+        # max number of directed outgoing arcs
+        n_max_dir_out = sum(
+            sum(
+                1
+                for j in m.set_J[(g, l, l_line)]
+                if j not in m.set_J_und[(g, l, l_line)]
+            )  # directed
+            for l_line in m.set_L[g]
+            if l_line != l
+            # if l_line not in m.set_L_exp[g] # cannot be an export: why?
+            if (g, l, l_line) in m.set_J
+        )
+        
+        # check the maximum number of incoming arcs
+        if n_max_dir_out <= 1:
+            # there can only be one outgoing arc at most: redundant constraint
+            # TODO: consider this condition when defining the set
+            return pyo.Constraint.Skip
+        else:  # more than one outgoing arc is possible
+
+            # number of (new) incoming directed arcs in a group
+            b_max_out_gl = 0
+
+            # the big m
+            M_gl = n_max_dir_out - 1  # has to be positive since n_max_dir_out > 1
+            # TODO: put parenthesis to avoid funny results
+            temp_constr = (
+                sum(
+                    # *********************************************************
+                    # interfaced groups
+                    sum(
+                        sum(
+                            1
+                            for j in m.set_J_col[(g, l, l_circ)]  # part of group
+                            if j not in m.set_J_und[(g, l, l_circ)]  # directed
+                            if (g, l, l_circ, j) in m.set_GLLJ_col_t[t]
+                        )
+                        * m.var_xi_arc_inv_t[t]  # in t
+                        for t in m.set_T_int
+                    )
+                    +
+                    # *********************************************************
+                    # optional non-interfaced groups
+                    sum(
+                        sum(
+                            sum(
+                                1
+                                for j in m.set_J_col[(g, l, l_circ)]  # part of group
+                                if j not in m.set_J_und[(g, l, l_circ)]  # directed
+                                if (g, l, l_circ, j) in m.set_GLLJ_col_t[t]
+                            )
+                            * m.var_delta_arc_inv_th[(t, h)]
+                            for h in m.set_H_t[t]
+                        )
+                        for t in m.set_T
+                        if t not in m.set_T_mdt  # optional
+                        if t not in m.set_T_int  # not interfaced
+                    )
+                    +
+                    # *********************************************************
+                    # interfaced arcs
+                    (sum(
+                        m.var_xi_arc_inv_gllj[(g, l, l_circ, j_circ)]
+                        for j_circ in m.set_J[(g, l, l_circ)]
+                        if j_circ not in m.set_J_und[(g, l, l_circ)]  # directed
+                        if j_circ in m.set_J_int[(g, l, l_circ)]  # interfaced
+                        if j_circ not in m.set_J_col[(g, l, l_circ)]  # individual
+                    )
+                    if (g, l, l_circ) in m.set_J
+                    else 0) +
+                    # *********************************************************
+                    # optional non-interfaced arcs
+                    (sum(
+                        sum(
+                            m.var_delta_arc_inv_glljh[(g, l, l_circ, j_dot, h_dot)]
+                            for h_dot in m.set_H_gllj[(g, l, l_circ, j_dot)]
+                        )
+                        for j_dot in m.set_J[(g, l, l_circ)]
+                        if j_dot not in m.set_J_und[(g, l, l_circ)]  # directed
+                        if j_dot not in m.set_J_int[(g, l, l_circ)]  # not interfaced
+                        if j_dot not in m.set_J_col[(g, l, l_circ)]  # individual
+                        if j_dot not in m.set_J_pre[(g, l, l_circ)]  # new
+                        if j_dot not in m.set_J_mdt[(g, l, l_circ)]  # optional
+                    )
+                    if (g, l, l_circ) in m.set_J
+                    else 0) +
+                    # *********************************************************
+                    # preexisting directed arcs
+                    (sum(
+                        1
+                        for j_pre_dir in m.set_J_pre[(g, l, l_circ)]  # preexisting
+                        if j_pre_dir not in m.set_J_und[(g, l, l_circ)]  # directed
+                    )
+                    if (g, l, l_circ) in m.set_J_pre
+                    else 0) +
+                    # *********************************************************
+                    # mandatory directed arcs
+                    (sum(
+                        1
+                        for j_mdt_dir in m.set_J_mdt[(g, l, l_circ)]
+                        if j_mdt_dir not in m.set_J_und[(g, l, l_circ)]  # directed
+                    )
+                    if (g, l, l_circ) in m.set_J_mdt
+                    else 0)
+                    # *********************************************************
+                    for l_circ in m.set_L[g]
+                    if l_circ not in m.set_L_imp[g]
+                    if l_circ != l
+                )
+                <= 1  # + 
+                # TODO: what is below has copy&pasted, must be completely revised
+                # M_gl*sum(
+                #     # *********************************************************
+                #     # outgoing arcs in interfaced groups, nominal direction
+                #     sum(sum(1
+                #             for j in m.set_J_col[(g,l,l_diamond)]
+                #             #if j in m.set_J_int[(g,l,l_diamond)]
+                #             if (g,l,l_diamond,j) in m.set_GLLJ_col_t[t]
+                #             )*m.var_xi_arc_inv_t[t]
+                #         for t in m.set_T_int
+                #         ) if (g,l,l_diamond) in m.set_J_col else 0
+                #     +
+                #     # outgoing arcs in interfaced groups, reverse direction
+                #     sum(sum(1
+                #             for j in m.set_J_col[(g,l_diamond,l)]
+                #             #if j in m.set_J_int[(g,l_diamond,l)]
+                #             if j in m.set_J_und[(g,l_diamond,l)]
+                #             if (g,l_diamond,l,j) in m.set_GLLJ_col_t[t]
+                #             )*m.var_xi_arc_inv_t[t]
+                #         for t in m.set_T_int
+                #         ) if (g,l_diamond,l) in m.set_J_col else 0
+                #     +
+                #     # *********************************************************
+                #     # TODO: outgoing arcs in non-interfaced optional groups, nominal
+                #     sum(sum(1
+                #             for j in m.set_J_col[(g,l,l_diamond)]
+                #             #if j in m.set_J_int[(g,l,l_diamond)]
+                #             if (g,l,l_diamond,j) in m.set_GLLJ_col_t[t]
+                #             )*sum(
+                #                 m.var_delta_arc_inv_th[(t,h)]
+                #                 for h in m.set_H_t[t]
+                #                 )
+                #         for t in m.set_T
+                #         if t not in m.set_T_mdt
+                #         if t not in m.set_T_int
+                #         ) if (g,l,l_diamond) in m.set_J_col else 0
+                #     +
+                #     # TODO: outgoing arcs in non-interfaced optional groups, reverse
+                #     sum(sum(1
+                #             for j in m.set_J_col[(g,l_diamond,l)]
+                #             #if j in m.set_J_int[(g,l_diamond,l)]
+                #             if j in m.set_J_und[(g,l_diamond,l)]
+                #             if (g,l_diamond,l,j) in m.set_GLLJ_col_t[t]
+                #             )*sum(
+                #                 m.var_delta_arc_inv_th[(t,h)]
+                #                 for h in m.set_H_t[t]
+                #                 )
+                #         for t in m.set_T
+                #         if t not in m.set_T_mdt
+                #         if t not in m.set_T_int
+                #         ) if (g,l_diamond,l) in m.set_J_col else 0
+                #     +
+                #     # *********************************************************
+                #     # interfaced individual outgoing arcs, nominal direction
+                #     sum(m.var_xi_arc_inv_gllj[(g,l,l_diamond,j)]
+                #         for j in m.set_J_int[(g,l,l_diamond)] # interfaced
+                #         if j not in m.set_J_col[(g,l,l_diamond)] # individual
+                #         ) if (g,l,l_diamond) in m.set_J_int else 0
+                #     +
+                #     # *********************************************************
+                #     # interfaced individual undirected arcs, reverse direction
+                #     sum(m.var_xi_arc_inv_gllj[(g,l,l_diamond,j)]
+                #         for j in m.set_J_und[(g,l_diamond,l)] # undirected
+                #         if j in m.set_J_int[(g,l_diamond,l)] # interfaced
+                #         if j not in m.set_J_col[(g,l_diamond,l)] # individual
+                #         ) if (g,l_diamond,l) in m.set_J_und else 0
+                #     +
+                #     # *********************************************************
+                #     # outgoing non-interfaced individual optional arcs
+                #     sum(
+                #         sum(m.var_delta_arc_inv_glljh[(g,l,l_diamond,j,h)]
+                #             for h in m.set_H_gllj[(g,l,l_diamond,j)])
+                #         for j in m.set_J[(g,l,l_diamond)]
+                #         if j not in m.set_J_col[(g,l,l_diamond)] # individual
+                #         if j not in m.set_J_mdt[(g,l,l_diamond)] # optional
+                #         if j not in m.set_J_int[(g,l,l_diamond)] # interfaced
+                #         ) if (g,l,l_diamond) in m.set_J else 0
+                #     +
+                #     # *********************************************************
+                #     # individual non-interfaced undirected arcs, reverse dir.
+                #     sum(
+                #         sum(m.var_delta_arc_inv_glljh[(g,l_diamond,l,j,h)]
+                #             for h in m.set_H_gllj[(g,l_diamond,l,j)])
+                #         for j in m.set_J_und[(g,l_diamond,l)] # undirected
+                #         if j not in m.set_J_col[(g,l_diamond,l)] # individual
+                #         if j not in m.set_J_mdt[(g,l_diamond,l)] # optional
+                #         if j not in m.set_J_int[(g,l_diamond,l)] # interfaced
+                #         ) if (g,l_diamond,l) in m.set_J_und else 0
+                #     +
+                #     # *********************************************************
+                #     # preselected outgonig arcs, nominal direction
+                #     len(m.set_J_pre[(g,l,l_diamond)]
+                #         ) if (g,l,l_diamond) in m.set_J_pre else 0
+                #     +
+                #     # *********************************************************
+                #     # mandatory outgoing arcs, nominal direction
+                #     len(m.set_J_mdt[(g,l,l_diamond)]
+                #         ) if (g,l,l_diamond) in m.set_J_mdt else 0
+                #     +
+                #     # *********************************************************
+                #     # undirected preselected arcs, reverse direction
+                #     sum(1
+                #         for j in m.set_J_pre[(g,l_diamond,l)]
+                #         if j in m.set_J_und[(g,l_diamond,l)]
+                #         ) if (g,l_diamond,l) in m.set_J_pre else 0
+                #     +
+                #     # *********************************************************
+                #     # undirected mandatory arcs, reverse direction
+                #     sum(1
+                #         for j in m.set_J_mdt[(g,l_diamond,l)]
+                #         if j in m.set_J_und[(g,l_diamond,l)]
+                #         ) if (g,l_diamond,l) in m.set_J_mdt else 0
+                #     # *********************************************************
+                #     for l_diamond in m.set_L[g]
+                #     if l_diamond not in m.set_L_imp[g]
+                #     if l_diamond != l
+                #     )
+            )
+            if type(temp_constr) == bool:
+                # trivial outcome
+                return pyo.Constraint.Feasible if temp_constr else pyo.Constraint.Infeasible
+            else:
+                # constraint is relevant
+                return temp_constr
+
+    model.constr_max_outgoing_directed_arcs = pyo.Constraint(
+        model.set_GL, rule=rule_constr_max_outgoing_directed_arcs
+    )
+
+    # *************************************************************************
+    # *************************************************************************
+
+    return model
+
+    # *************************************************************************
+    # *************************************************************************
+
+# *****************************************************************************
+# *****************************************************************************
diff --git a/src/topupopt/problems/esipp/blocks/prices.py b/src/topupopt/problems/esipp/blocks/prices.py
new file mode 100644
index 0000000..e648834
--- /dev/null
+++ b/src/topupopt/problems/esipp/blocks/prices.py
@@ -0,0 +1,205 @@
+# imports
+
+import pyomo.environ as pyo
+# *****************************************************************************
+# *****************************************************************************
+
+def add_prices_block(
+    model: pyo.AbstractModel,
+    enable_default_values: bool = True,
+    enable_validation: bool = True,
+    enable_initialisation: bool = True,
+):  
+    
+    # *************************************************************************
+    # *************************************************************************
+
+    # sparse index sets
+
+    # *************************************************************************
+
+    # set of price segments
+
+    model.set_S = pyo.Set(model.set_GL_exp_imp, model.set_QPK)
+
+    # set of GLQKS tuples
+
+    def init_set_GLQPKS(m):
+        return (
+            (g, l, q, p, k, s)
+            # for (g,l) in m.set_GL_exp_imp
+            # for (q,k) in m.set_QK
+            for (g, l, q, p, k) in m.set_S
+            for s in m.set_S[(g, l, q, p, k)]
+        )
+
+    model.set_GLQPKS = pyo.Set(
+        dimen=6, initialize=(init_set_GLQPKS if enable_initialisation else None)
+    )
+
+    def init_set_GLQPKS_exp(m):
+        return (
+            glqpks for glqpks in m.set_GLQPKS if glqpks[1] in m.set_L_exp[glqpks[0]]
+        )
+
+    model.set_GLQPKS_exp = pyo.Set(
+        dimen=6, initialize=(init_set_GLQPKS_exp if enable_initialisation else None)
+    )
+
+    def init_set_GLQPKS_imp(m):
+        return (
+            glqpks for glqpks in m.set_GLQPKS if glqpks[1] in m.set_L_imp[glqpks[0]]
+        )
+
+    model.set_GLQPKS_imp = pyo.Set(
+        dimen=6, initialize=(init_set_GLQPKS_imp if enable_initialisation else None)
+    )
+
+    # *************************************************************************
+    # *************************************************************************
+
+    # parameters
+
+    # *************************************************************************
+    # *************************************************************************
+
+    # objective function
+
+    # *************************************************************************
+
+    # resource prices
+
+    model.param_p_glqpks = pyo.Param(model.set_GLQPKS, within=pyo.NonNegativeReals)
+
+    # maximum resource volumes for each prices
+
+    model.param_v_max_glqpks = pyo.Param(model.set_GLQPKS, within=pyo.NonNegativeReals)
+
+    # *************************************************************************
+    # *************************************************************************
+
+    # variables
+
+    # *************************************************************************
+    # *************************************************************************
+
+    # objective function
+
+    # capex
+
+    # exported flow revenue
+
+    model.var_efr_glqpk = pyo.Var(
+        model.set_GL_exp, model.set_QPK, within=pyo.NonNegativeReals
+    )
+
+    # imported flow cost
+
+    model.var_ifc_glqpk = pyo.Var(
+        model.set_GL_imp, model.set_QPK, within=pyo.NonNegativeReals
+    )
+
+    # exported flow
+
+    # TODO: validate the bounds by ensuring inf. cap. only exists in last segm.
+
+    def bounds_var_ef_glqpks(m, g, l, q, p, k, s):
+        if (g, l, q, p, k, s) in m.param_v_max_glqpks:
+            # predefined finite capacity
+            return (0, m.param_v_max_glqpks[(g, l, q, p, k, s)])
+        else:
+            # infinite capacity
+            return (0, None)
+
+    model.var_ef_glqpks = pyo.Var(
+        model.set_GLQPKS_exp, within=pyo.NonNegativeReals, bounds=bounds_var_ef_glqpks
+    )
+
+    # imported flow
+
+    def bounds_var_if_glqpks(m, g, l, q, p, k, s):
+        if (g, l, q, p, k, s) in m.param_v_max_glqpks:
+            # predefined finite capacity
+            return (0, m.param_v_max_glqpks[(g, l, q, p, k, s)])
+        else:
+            # infinite capacity
+            return (0, None)
+
+    model.var_if_glqpks = pyo.Var(
+        model.set_GLQPKS_imp, within=pyo.NonNegativeReals, bounds=bounds_var_if_glqpks
+    )
+
+    # *************************************************************************
+    # *************************************************************************
+
+    # exported flow revenue
+    def rule_constr_exp_flow_revenue(m, g, l, q, p, k):
+        return (
+            sum(
+                m.var_ef_glqpks[(g, l, q, p, k, s)]
+                * m.param_p_glqpks[(g, l, q, p, k, s)]
+                for s in m.set_S[(g, l, q, p, k)]
+            )
+            == m.var_efr_glqpk[(g, l, q, p, k)]
+        )
+
+    model.constr_exp_flow_revenue = pyo.Constraint(
+        model.set_GL_exp, model.set_QPK, rule=rule_constr_exp_flow_revenue
+    )
+
+    # imported flow cost
+    def rule_constr_imp_flow_cost(m, g, l, q, p, k):
+        return (
+            sum(
+                m.var_if_glqpks[(g, l, q, p, k, s)]
+                * m.param_p_glqpks[(g, l, q, p, k, s)]
+                for s in m.set_S[(g, l, q, p, k)]
+            )
+            == m.var_ifc_glqpk[(g, l, q, p, k)]
+        )
+
+    model.constr_imp_flow_cost = pyo.Constraint(
+        model.set_GL_imp, model.set_QPK, rule=rule_constr_imp_flow_cost
+    )
+
+    # exported flows
+    def rule_constr_exp_flows(m, g, l, q, p, k):
+        return sum(
+            m.var_v_glljqk[(g, l_star, l, j, q, k)]
+            * m.param_eta_glljqk[(g, l_star, l, j, q, k)]
+            for l_star in m.set_L[g]
+            if l_star not in m.set_L_exp[g]
+            for j in m.set_J[(g, l_star, l)]  # only directed arcs
+        ) == sum(m.var_ef_glqpks[(g, l, q, p, k, s)] for s in m.set_S[(g, l, q, p, k)])
+
+    model.constr_exp_flows = pyo.Constraint(
+        model.set_GL_exp, model.set_QPK, rule=rule_constr_exp_flows
+    )
+
+    # imported flows
+    def rule_constr_imp_flows(m, g, l, q, p, k):
+        return sum(
+            m.var_v_glljqk[(g, l, l_star, j, q, k)]
+            for l_star in m.set_L[g]
+            if l_star not in m.set_L_imp[g]
+            for j in m.set_J[(g, l, l_star)]  # only directed arcs
+        ) == sum(m.var_if_glqpks[(g, l, q, p, k, s)] for s in m.set_S[(g, l, q, p, k)])
+
+    model.constr_imp_flows = pyo.Constraint(
+        model.set_GL_imp, model.set_QPK, rule=rule_constr_imp_flows
+    )
+    
+    # *************************************************************************
+    # *************************************************************************
+
+    return model
+
+    # *************************************************************************
+    # *************************************************************************
+
+
+# *****************************************************************************
+# *****************************************************************************
+# *****************************************************************************
+# *****************************************************************************
+# *****************************************************************************
diff --git a/src/topupopt/problems/esipp/model.py b/src/topupopt/problems/esipp/model.py
index 412e5d9..566e228 100644
--- a/src/topupopt/problems/esipp/model.py
+++ b/src/topupopt/problems/esipp/model.py
@@ -2,7 +2,7 @@
 
 import pyomo.environ as pyo
 
-from math import isfinite, inf
+from .blocks.networks import add_network_restrictions
 
 # *****************************************************************************
 # *****************************************************************************
@@ -85,17 +85,6 @@ def create_model(
 
     model.set_L_exp = pyo.Set(model.set_G, within=model.set_L)
 
-    # set of nodes on network g incompatible with having more than one incoming
-    # arc unless there are outgoing arcs too
-
-    model.set_L_max_in_g = pyo.Set(
-        model.set_G, within=model.set_L
-    )  # should inherently exclude import nodes
-
-    model.set_L_max_out_g = pyo.Set(
-        model.set_G, within=model.set_L
-    )  # should inherently exclude export nodes
-
     # *************************************************************************
     # *************************************************************************
 
@@ -1621,26 +1610,6 @@ def create_model(
         model.set_GL_not_exp_imp, model.set_QK, within=pyo.Reals, default=0
     )
 
-    # maximum number of arcs per node pair
-
-    model.param_max_number_parallel_arcs = pyo.Param(
-        model.set_GLL,
-        # within=pyo.PositiveIntegers,
-        within=pyo.PositiveReals,
-        default=inf,
-    )
-
-    def init_set_GLL_arc_max(m):
-        return (
-            (g, l1, l2)
-            for (g, l1, l2) in m.param_max_number_parallel_arcs
-            if isfinite(m.param_max_number_parallel_arcs[(g, l1, l2)])
-        )
-
-    model.set_GLL_arc_max = pyo.Set(
-        dimen=3, within=model.set_GLL, initialize=init_set_GLL_arc_max
-    )
-
     # effect of system inputs on specific network and node pairs
 
     model.param_a_nw_glimqk = pyo.Param(
@@ -2479,806 +2448,9 @@ def create_model(
     )
 
     # *************************************************************************
-
-    # limit number of directed arcs per direction
-
-    def rule_constr_limited_parallel_arcs_per_direction(m, g, l1, l2):
-        # cases:
-        # 1) the number of options is lower than or equal to the limit (skip)
-        # 2) the number of preexisting and new mandatory arcs exceeds
-        # the limit (infeasible: pyo.Constraint.Infeasible)
-        # 3) all other cases (constraint)
-
-        # number of preexisting arcs going from l1 to l2
-
-        number_arcs_pre_nom = (
-            len(m.set_J_pre[(g, l1, l2)]) if (g, l1, l2) in m.set_J_pre else 0
-        )
-
-        number_arcs_pre_rev = (
-            sum(1 for j in m.set_J_pre[(g, l2, l1)] if j in m.set_J_und[(g, l2, l1)])
-            if (g, l2, l1) in m.set_J_pre
-            else 0
-        )
-
-        # number of mandatory arcs going from l1 to l2
-
-        number_arcs_mdt_nom = (
-            len(m.set_J_mdt[(g, l1, l2)]) if (g, l1, l2) in m.set_J_mdt else 0
-        )
-
-        number_arcs_mdt_rev = (
-            sum(1 for j in m.set_J_mdt[(g, l2, l1)] if j in m.set_J_und[(g, l2, l1)])
-            if (g, l2, l1) in m.set_J_mdt
-            else 0
-        )
-
-        # number of optional arcs going from l1 to l2
-
-        number_arcs_opt_nom = (
-            sum(
-                1
-                for j in m.set_J[(g, l1, l2)]
-                if j not in m.set_J_pre[(g, l1, l2)]
-                if j not in m.set_J_mdt[(g, l1, l2)]
-            )
-            if (g, l1, l2) in m.set_J
-            else 0
-        )
-
-        number_arcs_opt_rev = (
-            sum(
-                1
-                for j in m.set_J[(g, l2, l1)]
-                if j not in m.set_J_pre[(g, l2, l1)]
-                if j not in m.set_J_mdt[(g, l2, l1)]
-                if j in m.set_J_und[(g, l2, l1)]
-            )
-            if (g, l2, l1) in m.set_J
-            else 0
-        )
-
-        # build the constraints
-
-        if (
-            number_arcs_mdt_nom
-            + number_arcs_mdt_rev
-            + number_arcs_pre_nom
-            + number_arcs_pre_rev
-            > m.param_max_number_parallel_arcs[(g, l1, l2)]
-        ):
-            # the number of unavoidable arcs already exceeds the limit
-
-            return pyo.Constraint.Infeasible
-
-        elif (
-            number_arcs_opt_nom
-            + number_arcs_opt_rev
-            + number_arcs_mdt_nom
-            + number_arcs_mdt_rev
-            + number_arcs_pre_nom
-            + number_arcs_pre_rev
-            > m.param_max_number_parallel_arcs[(g, l1, l2)]
-        ):
-            # the number of potential arcs exceeds the limit: cannot be skipped
-
-            return (
-                # preexisting arcs
-                number_arcs_pre_nom + number_arcs_pre_rev +
-                # mandatory arcs
-                number_arcs_mdt_nom + number_arcs_mdt_rev +
-                # arcs within an (optional) group that uses interfaces
-                sum(
-                    (
-                        sum(
-                            1
-                            for j in m.set_J_col[(g, l1, l2)]
-                            if (g, l1, l2, j) in m.set_GLLJ_col_t[t]
-                        )
-                        if (g, l1, l2) in m.set_J_col
-                        else 0
-                        + sum(
-                            1
-                            for j in m.set_J_col[(g, l2, l1)]
-                            if j in m.set_J_und[(g, l2, l1)]
-                            if (g, l2, l1, j) in m.set_GLLJ_col_t[t]
-                        )
-                        if ((g, l2, l1) in m.set_J_col and (g, l2, l1) in m.set_J_und)
-                        else 0
-                    )
-                    * m.var_xi_arc_inv_t[t]
-                    for t in m.set_T_int
-                )
-                +
-                # arcs within an (optional) group that does not use interfaces
-                sum(
-                    (
-                        sum(
-                            1
-                            for j in m.set_J_col[(g, l1, l2)]
-                            if (g, l1, l2, j) in m.set_GLLJ_col_t[t]
-                        )
-                        if (g, l1, l2) in m.set_J_col
-                        else 0
-                        + sum(
-                            1
-                            for j in m.set_J_col[(g, l2, l1)]
-                            if j in m.set_J_und[(g, l2, l1)]
-                            if (g, l2, l1, j) in m.set_GLLJ_col_t[t]
-                        )
-                        if ((g, l2, l1) in m.set_J_col and (g, l2, l1) in m.set_J_und)
-                        else 0
-                    )
-                    * sum(m.var_delta_arc_inv_th[(t, h)] for h in m.set_H_t[t])
-                    for t in m.set_T  # new
-                    if t not in m.set_T_mdt  # optional
-                    if t not in m.set_T_int  # not interfaced
-                )
-                +
-                # optional individual arcs using interfaces, nominal direction
-                sum(
-                    m.var_xi_arc_inv_gllj[(g, l1, l2, j)]
-                    for j in m.set_J_int[(g, l1, l2)]  # interfaced
-                    if j not in m.set_J_col[(g, l1, l2)]  # individual
-                )
-                if (g, l1, l2) in m.set_J_int
-                else 0 +
-                # optional individual arcs using interfaces, reverse direction
-                sum(
-                    m.var_xi_arc_inv_gllj[(g, l2, l1, j)]
-                    for j in m.set_J_int[(g, l2, l1)]  # interfaced
-                    if j in m.set_J_und[(g, l2, l1)]  # undirected
-                    if j not in m.set_J_col[(g, l1, l2)]  # individual
-                )
-                if ((g, l2, l1) in m.set_J_int and (g, l2, l1) in m.set_J_und)
-                else 0 +
-                # optional individual arcs not using interfaces, nominal dir.
-                sum(
-                    sum(
-                        m.var_delta_arc_inv_glljh[(g, l1, l2, j, h)]
-                        for h in m.set_H_gllj[(g, l1, l2, j)]
-                    )
-                    for j in m.set_J[(g, l1, l2)]
-                    if j not in m.set_J_pre[(g, l1, l2)]  # not preexisting
-                    if j not in m.set_J_mdt[(g, l1, l2)]  # not mandatory
-                    if j not in m.set_J_int[(g, l1, l2)]  # not interfaced
-                    if j not in m.set_J_col[(g, l1, l2)]  # individual
-                )
-                if (g, l1, l2) in m.set_J
-                else 0 +
-                # optional individual arcs not using interfaces, reverse dir.
-                sum(
-                    sum(
-                        m.var_delta_arc_inv_glljh[(g, l2, l1, j, h)]
-                        for h in m.set_H_gllj[(g, l2, l1, j)]
-                    )
-                    for j in m.set_J_opt[(g, l2, l1)]
-                    if j in m.set_J_und[(g, l2, l1)]
-                    if j not in m.set_J_pre[(g, l2, l1)]  # not preexisting
-                    if j not in m.set_J_mdt[(g, l2, l1)]  # not mandatory
-                    if j not in m.set_J_int[(g, l2, l1)]  # not interfaced
-                    if j not in m.set_J_col[(g, l2, l1)]  # individual
-                )
-                if (g, l2, l1) in m.set_J
-                else 0 <= m.param_max_number_parallel_arcs[(g, l1, l2)]
-            )
-
-        else:  # the number of options is lower than or equal to the limit: skip
-            return pyo.Constraint.Skip
-
-    model.constr_limited_parallel_arcs_per_direction = pyo.Constraint(
-        model.set_GLL_arc_max, rule=rule_constr_limited_parallel_arcs_per_direction
-    )
-
-    # *************************************************************************
-
-    # there can only one incoming arc at most, if there are no outgoing arcs
-
-    def rule_constr_max_incoming_directed_arcs(m, g, l):
-        # check if the node is not among those subject to a limited number of incoming arcs
-        if l not in m.set_L_max_in_g[g]:
-            # it is not, skip this constraint
-            return pyo.Constraint.Skip
-
-        # max number of directed incoming arcs
-        n_max_dir_in = sum(
-            sum(
-                1
-                for j in m.set_J[(g, l_line, l)]
-                if j not in m.set_J_und[(g, l_line, l)]
-            )  # directed
-            for l_line in m.set_L[g] # for every node
-            if l_line != l # cannot be the same node
-            # if l_line not in m.set_L_imp[g] # why?
-            if (g, l_line, l) in m.set_J
-        )
-        
-        # check the maximum number of incoming arcs
-        if n_max_dir_in <= 1:
-            # there can only be one incoming arc at most: redundant constraint
-            return pyo.Constraint.Skip
-        else:  # more than one incoming arc is possible
-        
-            # number of (new) incoming directed arcs in a group
-            b_max_in_gl = 0
-
-            # the big m
-
-            M_gl = n_max_dir_in - 1  # has to be positive since n_max_dir_in > 1
-            # TODO: put parenthesis to avoid funny results
-            temp_constr = (
-                sum(
-                    # *********************************************************
-                    # interfaced groups
-                    sum(
-                        sum(
-                            1
-                            for j in m.set_J_col[(g, l_circ, l)]  # part of group
-                            if j not in m.set_J_und[(g, l_circ, l)]  # directed
-                            if (g, l_circ, l, j) in m.set_GLLJ_col_t[t]
-                        )
-                        * m.var_xi_arc_inv_t[t]  # in t
-                        for t in m.set_T_int
-                    )
-                    +
-                    # *********************************************************
-                    # optional non-interfaced groups
-                    sum(
-                        sum(
-                            sum(
-                                1
-                                for j in m.set_J_col[(g, l_circ, l)]  # part of group
-                                if j not in m.set_J_und[(g, l_circ, l)]  # directed
-                                if (g, l_circ, l, j) in m.set_GLLJ_col_t[t]
-                            )
-                            * m.var_delta_arc_inv_th[(t, h)]
-                            for h in m.set_H_t[t]
-                        )
-                        for t in m.set_T
-                        if t not in m.set_T_mdt  # optional
-                        if t not in m.set_T_int  # not interfaced
-                    )
-                    +
-                    # *********************************************************
-                    # interfaced arcs
-                    (sum(
-                        m.var_xi_arc_inv_gllj[(g, l_circ, l, j_circ)]
-                        for j_circ in m.set_J[(g, l_circ, l)]
-                        if j_circ not in m.set_J_und[(g, l_circ, l)]  # directed
-                        if j_circ in m.set_J_int[(g, l_circ, l)]  # interfaced
-                        if j_circ not in m.set_J_col[(g, l_circ, l)]  # individual
-                    )
-                    if (g, l_circ, l) in m.set_J
-                    else 0) +
-                    # *********************************************************
-                    # optional non-interfaced arcs
-                    (sum(
-                        sum(
-                            m.var_delta_arc_inv_glljh[(g, l_circ, l, j_dot, h_dot)]
-                            for h_dot in m.set_H_gllj[(g, l_circ, l, j_dot)]
-                        )
-                        for j_dot in m.set_J[(g, l_circ, l)]
-                        if j_dot not in m.set_J_und[(g, l_circ, l)]  # directed
-                        if j_dot not in m.set_J_int[(g, l_circ, l)]  # not interfaced
-                        if j_dot not in m.set_J_col[(g, l_circ, l)]  # individual
-                        if j_dot not in m.set_J_pre[(g, l_circ, l)]  # new
-                        if j_dot not in m.set_J_mdt[(g, l_circ, l)]  # optional
-                    )
-                    if (g, l_circ, l) in m.set_J
-                    else 0) +
-                    # *********************************************************
-                    # preexisting directed arcs
-                    (sum(
-                        1
-                        for j_pre_dir in m.set_J_pre[(g, l_circ, l)]  # preexisting
-                        if j_pre_dir not in m.set_J_und[(g, l_circ, l)]  # directed
-                    )
-                    if (g, l_circ, l) in m.set_J_pre
-                    else 0) +
-                    # *********************************************************
-                    # mandatory directed arcs
-                    (sum(
-                        1
-                        for j_mdt_dir in m.set_J_mdt[(g, l_circ, l)]
-                        if j_mdt_dir not in m.set_J_und[(g, l_circ, l)]  # directed
-                    )
-                    if (g, l_circ, l) in m.set_J_mdt
-                    else 0)
-                    # *********************************************************
-                    for l_circ in m.set_L[g]
-                    if l_circ not in m.set_L_exp[g]
-                    if l_circ != l
-                )
-                <= 1  # +
-                # M_gl*sum(
-                #     # *********************************************************
-                #     # outgoing arcs in interfaced groups, nominal direction
-                #     sum(sum(1
-                #             for j in m.set_J_col[(g,l,l_diamond)]
-                #             #if j in m.set_J_int[(g,l,l_diamond)]
-                #             if (g,l,l_diamond,j) in m.set_GLLJ_col_t[t]
-                #             )*m.var_xi_arc_inv_t[t]
-                #         for t in m.set_T_int
-                #         ) if (g,l,l_diamond) in m.set_J_col else 0
-                #     +
-                #     # outgoing arcs in interfaced groups, reverse direction
-                #     sum(sum(1
-                #             for j in m.set_J_col[(g,l_diamond,l)]
-                #             #if j in m.set_J_int[(g,l_diamond,l)]
-                #             if j in m.set_J_und[(g,l_diamond,l)]
-                #             if (g,l_diamond,l,j) in m.set_GLLJ_col_t[t]
-                #             )*m.var_xi_arc_inv_t[t]
-                #         for t in m.set_T_int
-                #         ) if (g,l_diamond,l) in m.set_J_col else 0
-                #     +
-                #     # *********************************************************
-                #     # TODO: outgoing arcs in non-interfaced optional groups, nominal
-                #     sum(sum(1
-                #             for j in m.set_J_col[(g,l,l_diamond)]
-                #             #if j in m.set_J_int[(g,l,l_diamond)]
-                #             if (g,l,l_diamond,j) in m.set_GLLJ_col_t[t]
-                #             )*sum(
-                #                 m.var_delta_arc_inv_th[(t,h)]
-                #                 for h in m.set_H_t[t]
-                #                 )
-                #         for t in m.set_T
-                #         if t not in m.set_T_mdt
-                #         if t not in m.set_T_int
-                #         ) if (g,l,l_diamond) in m.set_J_col else 0
-                #     +
-                #     # TODO: outgoing arcs in non-interfaced optional groups, reverse
-                #     sum(sum(1
-                #             for j in m.set_J_col[(g,l_diamond,l)]
-                #             #if j in m.set_J_int[(g,l_diamond,l)]
-                #             if j in m.set_J_und[(g,l_diamond,l)]
-                #             if (g,l_diamond,l,j) in m.set_GLLJ_col_t[t]
-                #             )*sum(
-                #                 m.var_delta_arc_inv_th[(t,h)]
-                #                 for h in m.set_H_t[t]
-                #                 )
-                #         for t in m.set_T
-                #         if t not in m.set_T_mdt
-                #         if t not in m.set_T_int
-                #         ) if (g,l_diamond,l) in m.set_J_col else 0
-                #     +
-                #     # *********************************************************
-                #     # interfaced individual outgoing arcs, nominal direction
-                #     sum(m.var_xi_arc_inv_gllj[(g,l,l_diamond,j)]
-                #         for j in m.set_J_int[(g,l,l_diamond)] # interfaced
-                #         if j not in m.set_J_col[(g,l,l_diamond)] # individual
-                #         ) if (g,l,l_diamond) in m.set_J_int else 0
-                #     +
-                #     # *********************************************************
-                #     # interfaced individual undirected arcs, reverse direction
-                #     sum(m.var_xi_arc_inv_gllj[(g,l,l_diamond,j)]
-                #         for j in m.set_J_und[(g,l_diamond,l)] # undirected
-                #         if j in m.set_J_int[(g,l_diamond,l)] # interfaced
-                #         if j not in m.set_J_col[(g,l_diamond,l)] # individual
-                #         ) if (g,l_diamond,l) in m.set_J_und else 0
-                #     +
-                #     # *********************************************************
-                #     # outgoing non-interfaced individual optional arcs
-                #     sum(
-                #         sum(m.var_delta_arc_inv_glljh[(g,l,l_diamond,j,h)]
-                #             for h in m.set_H_gllj[(g,l,l_diamond,j)])
-                #         for j in m.set_J[(g,l,l_diamond)]
-                #         if j not in m.set_J_col[(g,l,l_diamond)] # individual
-                #         if j not in m.set_J_mdt[(g,l,l_diamond)] # optional
-                #         if j not in m.set_J_int[(g,l,l_diamond)] # interfaced
-                #         ) if (g,l,l_diamond) in m.set_J else 0
-                #     +
-                #     # *********************************************************
-                #     # individual non-interfaced undirected arcs, reverse dir.
-                #     sum(
-                #         sum(m.var_delta_arc_inv_glljh[(g,l_diamond,l,j,h)]
-                #             for h in m.set_H_gllj[(g,l_diamond,l,j)])
-                #         for j in m.set_J_und[(g,l_diamond,l)] # undirected
-                #         if j not in m.set_J_col[(g,l_diamond,l)] # individual
-                #         if j not in m.set_J_mdt[(g,l_diamond,l)] # optional
-                #         if j not in m.set_J_int[(g,l_diamond,l)] # interfaced
-                #         ) if (g,l_diamond,l) in m.set_J_und else 0
-                #     +
-                #     # *********************************************************
-                #     # preselected outgonig arcs, nominal direction
-                #     len(m.set_J_pre[(g,l,l_diamond)]
-                #         ) if (g,l,l_diamond) in m.set_J_pre else 0
-                #     +
-                #     # *********************************************************
-                #     # mandatory outgoing arcs, nominal direction
-                #     len(m.set_J_mdt[(g,l,l_diamond)]
-                #         ) if (g,l,l_diamond) in m.set_J_mdt else 0
-                #     +
-                #     # *********************************************************
-                #     # undirected preselected arcs, reverse direction
-                #     sum(1
-                #         for j in m.set_J_pre[(g,l_diamond,l)]
-                #         if j in m.set_J_und[(g,l_diamond,l)]
-                #         ) if (g,l_diamond,l) in m.set_J_pre else 0
-                #     +
-                #     # *********************************************************
-                #     # undirected mandatory arcs, reverse direction
-                #     sum(1
-                #         for j in m.set_J_mdt[(g,l_diamond,l)]
-                #         if j in m.set_J_und[(g,l_diamond,l)]
-                #         ) if (g,l_diamond,l) in m.set_J_mdt else 0
-                #     # *********************************************************
-                #     for l_diamond in m.set_L[g]
-                #     if l_diamond not in m.set_L_imp[g]
-                #     if l_diamond != l
-                #     )
-            )
-            if type(temp_constr) == bool:
-                # trivial outcome
-                return pyo.Constraint.Feasible if temp_constr else pyo.Constraint.Infeasible
-            else:
-                # constraint is relevant
-                return temp_constr
-
-    model.constr_max_incoming_directed_arcs = pyo.Constraint(
-        model.set_GL, rule=rule_constr_max_incoming_directed_arcs
-    )
-
     # *************************************************************************
-
-    def rule_constr_max_outgoing_directed_arcs(m, g, l):
-        # check if the node is not among those subject to a limited number of outgoing arcs
-        if l not in m.set_L_max_out_g[g]:
-            # it is not, skip this constraint
-            return pyo.Constraint.Skip
-
-        # max number of directed outgoing arcs
-        n_max_dir_out = sum(
-            sum(
-                1
-                for j in m.set_J[(g, l, l_line)]
-                if j not in m.set_J_und[(g, l, l_line)]
-            )  # directed
-            for l_line in m.set_L[g]
-            if l_line != l
-            # if l_line not in m.set_L_exp[g] # cannot be an export: why?
-            if (g, l, l_line) in m.set_J
-        )
-        
-        # check the maximum number of incoming arcs
-        if n_max_dir_out <= 1:
-            # there can only be one outgoing arc at most: redundant constraint
-            # TODO: consider this condition when defining the set
-            return pyo.Constraint.Skip
-        else:  # more than one outgoing arc is possible
-
-            # number of (new) incoming directed arcs in a group
-            b_max_out_gl = 0
-
-            # the big m
-            M_gl = n_max_dir_out - 1  # has to be positive since n_max_dir_out > 1
-            # TODO: put parenthesis to avoid funny results
-            temp_constr = (
-                sum(
-                    # *********************************************************
-                    # interfaced groups
-                    sum(
-                        sum(
-                            1
-                            for j in m.set_J_col[(g, l, l_circ)]  # part of group
-                            if j not in m.set_J_und[(g, l, l_circ)]  # directed
-                            if (g, l, l_circ, j) in m.set_GLLJ_col_t[t]
-                        )
-                        * m.var_xi_arc_inv_t[t]  # in t
-                        for t in m.set_T_int
-                    )
-                    +
-                    # *********************************************************
-                    # optional non-interfaced groups
-                    sum(
-                        sum(
-                            sum(
-                                1
-                                for j in m.set_J_col[(g, l, l_circ)]  # part of group
-                                if j not in m.set_J_und[(g, l, l_circ)]  # directed
-                                if (g, l, l_circ, j) in m.set_GLLJ_col_t[t]
-                            )
-                            * m.var_delta_arc_inv_th[(t, h)]
-                            for h in m.set_H_t[t]
-                        )
-                        for t in m.set_T
-                        if t not in m.set_T_mdt  # optional
-                        if t not in m.set_T_int  # not interfaced
-                    )
-                    +
-                    # *********************************************************
-                    # interfaced arcs
-                    (sum(
-                        m.var_xi_arc_inv_gllj[(g, l, l_circ, j_circ)]
-                        for j_circ in m.set_J[(g, l, l_circ)]
-                        if j_circ not in m.set_J_und[(g, l, l_circ)]  # directed
-                        if j_circ in m.set_J_int[(g, l, l_circ)]  # interfaced
-                        if j_circ not in m.set_J_col[(g, l, l_circ)]  # individual
-                    )
-                    if (g, l, l_circ) in m.set_J
-                    else 0) +
-                    # *********************************************************
-                    # optional non-interfaced arcs
-                    (sum(
-                        sum(
-                            m.var_delta_arc_inv_glljh[(g, l, l_circ, j_dot, h_dot)]
-                            for h_dot in m.set_H_gllj[(g, l, l_circ, j_dot)]
-                        )
-                        for j_dot in m.set_J[(g, l, l_circ)]
-                        if j_dot not in m.set_J_und[(g, l, l_circ)]  # directed
-                        if j_dot not in m.set_J_int[(g, l, l_circ)]  # not interfaced
-                        if j_dot not in m.set_J_col[(g, l, l_circ)]  # individual
-                        if j_dot not in m.set_J_pre[(g, l, l_circ)]  # new
-                        if j_dot not in m.set_J_mdt[(g, l, l_circ)]  # optional
-                    )
-                    if (g, l, l_circ) in m.set_J
-                    else 0) +
-                    # *********************************************************
-                    # preexisting directed arcs
-                    (sum(
-                        1
-                        for j_pre_dir in m.set_J_pre[(g, l, l_circ)]  # preexisting
-                        if j_pre_dir not in m.set_J_und[(g, l, l_circ)]  # directed
-                    )
-                    if (g, l, l_circ) in m.set_J_pre
-                    else 0) +
-                    # *********************************************************
-                    # mandatory directed arcs
-                    (sum(
-                        1
-                        for j_mdt_dir in m.set_J_mdt[(g, l, l_circ)]
-                        if j_mdt_dir not in m.set_J_und[(g, l, l_circ)]  # directed
-                    )
-                    if (g, l, l_circ) in m.set_J_mdt
-                    else 0)
-                    # *********************************************************
-                    for l_circ in m.set_L[g]
-                    if l_circ not in m.set_L_imp[g]
-                    if l_circ != l
-                )
-                <= 1  # + 
-                # TODO: what is below has copy&pasted, must be completely revised
-                # M_gl*sum(
-                #     # *********************************************************
-                #     # outgoing arcs in interfaced groups, nominal direction
-                #     sum(sum(1
-                #             for j in m.set_J_col[(g,l,l_diamond)]
-                #             #if j in m.set_J_int[(g,l,l_diamond)]
-                #             if (g,l,l_diamond,j) in m.set_GLLJ_col_t[t]
-                #             )*m.var_xi_arc_inv_t[t]
-                #         for t in m.set_T_int
-                #         ) if (g,l,l_diamond) in m.set_J_col else 0
-                #     +
-                #     # outgoing arcs in interfaced groups, reverse direction
-                #     sum(sum(1
-                #             for j in m.set_J_col[(g,l_diamond,l)]
-                #             #if j in m.set_J_int[(g,l_diamond,l)]
-                #             if j in m.set_J_und[(g,l_diamond,l)]
-                #             if (g,l_diamond,l,j) in m.set_GLLJ_col_t[t]
-                #             )*m.var_xi_arc_inv_t[t]
-                #         for t in m.set_T_int
-                #         ) if (g,l_diamond,l) in m.set_J_col else 0
-                #     +
-                #     # *********************************************************
-                #     # TODO: outgoing arcs in non-interfaced optional groups, nominal
-                #     sum(sum(1
-                #             for j in m.set_J_col[(g,l,l_diamond)]
-                #             #if j in m.set_J_int[(g,l,l_diamond)]
-                #             if (g,l,l_diamond,j) in m.set_GLLJ_col_t[t]
-                #             )*sum(
-                #                 m.var_delta_arc_inv_th[(t,h)]
-                #                 for h in m.set_H_t[t]
-                #                 )
-                #         for t in m.set_T
-                #         if t not in m.set_T_mdt
-                #         if t not in m.set_T_int
-                #         ) if (g,l,l_diamond) in m.set_J_col else 0
-                #     +
-                #     # TODO: outgoing arcs in non-interfaced optional groups, reverse
-                #     sum(sum(1
-                #             for j in m.set_J_col[(g,l_diamond,l)]
-                #             #if j in m.set_J_int[(g,l_diamond,l)]
-                #             if j in m.set_J_und[(g,l_diamond,l)]
-                #             if (g,l_diamond,l,j) in m.set_GLLJ_col_t[t]
-                #             )*sum(
-                #                 m.var_delta_arc_inv_th[(t,h)]
-                #                 for h in m.set_H_t[t]
-                #                 )
-                #         for t in m.set_T
-                #         if t not in m.set_T_mdt
-                #         if t not in m.set_T_int
-                #         ) if (g,l_diamond,l) in m.set_J_col else 0
-                #     +
-                #     # *********************************************************
-                #     # interfaced individual outgoing arcs, nominal direction
-                #     sum(m.var_xi_arc_inv_gllj[(g,l,l_diamond,j)]
-                #         for j in m.set_J_int[(g,l,l_diamond)] # interfaced
-                #         if j not in m.set_J_col[(g,l,l_diamond)] # individual
-                #         ) if (g,l,l_diamond) in m.set_J_int else 0
-                #     +
-                #     # *********************************************************
-                #     # interfaced individual undirected arcs, reverse direction
-                #     sum(m.var_xi_arc_inv_gllj[(g,l,l_diamond,j)]
-                #         for j in m.set_J_und[(g,l_diamond,l)] # undirected
-                #         if j in m.set_J_int[(g,l_diamond,l)] # interfaced
-                #         if j not in m.set_J_col[(g,l_diamond,l)] # individual
-                #         ) if (g,l_diamond,l) in m.set_J_und else 0
-                #     +
-                #     # *********************************************************
-                #     # outgoing non-interfaced individual optional arcs
-                #     sum(
-                #         sum(m.var_delta_arc_inv_glljh[(g,l,l_diamond,j,h)]
-                #             for h in m.set_H_gllj[(g,l,l_diamond,j)])
-                #         for j in m.set_J[(g,l,l_diamond)]
-                #         if j not in m.set_J_col[(g,l,l_diamond)] # individual
-                #         if j not in m.set_J_mdt[(g,l,l_diamond)] # optional
-                #         if j not in m.set_J_int[(g,l,l_diamond)] # interfaced
-                #         ) if (g,l,l_diamond) in m.set_J else 0
-                #     +
-                #     # *********************************************************
-                #     # individual non-interfaced undirected arcs, reverse dir.
-                #     sum(
-                #         sum(m.var_delta_arc_inv_glljh[(g,l_diamond,l,j,h)]
-                #             for h in m.set_H_gllj[(g,l_diamond,l,j)])
-                #         for j in m.set_J_und[(g,l_diamond,l)] # undirected
-                #         if j not in m.set_J_col[(g,l_diamond,l)] # individual
-                #         if j not in m.set_J_mdt[(g,l_diamond,l)] # optional
-                #         if j not in m.set_J_int[(g,l_diamond,l)] # interfaced
-                #         ) if (g,l_diamond,l) in m.set_J_und else 0
-                #     +
-                #     # *********************************************************
-                #     # preselected outgonig arcs, nominal direction
-                #     len(m.set_J_pre[(g,l,l_diamond)]
-                #         ) if (g,l,l_diamond) in m.set_J_pre else 0
-                #     +
-                #     # *********************************************************
-                #     # mandatory outgoing arcs, nominal direction
-                #     len(m.set_J_mdt[(g,l,l_diamond)]
-                #         ) if (g,l,l_diamond) in m.set_J_mdt else 0
-                #     +
-                #     # *********************************************************
-                #     # undirected preselected arcs, reverse direction
-                #     sum(1
-                #         for j in m.set_J_pre[(g,l_diamond,l)]
-                #         if j in m.set_J_und[(g,l_diamond,l)]
-                #         ) if (g,l_diamond,l) in m.set_J_pre else 0
-                #     +
-                #     # *********************************************************
-                #     # undirected mandatory arcs, reverse direction
-                #     sum(1
-                #         for j in m.set_J_mdt[(g,l_diamond,l)]
-                #         if j in m.set_J_und[(g,l_diamond,l)]
-                #         ) if (g,l_diamond,l) in m.set_J_mdt else 0
-                #     # *********************************************************
-                #     for l_diamond in m.set_L[g]
-                #     if l_diamond not in m.set_L_imp[g]
-                #     if l_diamond != l
-                #     )
-            )
-            if type(temp_constr) == bool:
-                # trivial outcome
-                return pyo.Constraint.Feasible if temp_constr else pyo.Constraint.Infeasible
-            else:
-                # constraint is relevant
-                return temp_constr
-
-    model.constr_max_outgoing_directed_arcs = pyo.Constraint(
-        model.set_GL, rule=rule_constr_max_outgoing_directed_arcs
-    )
-
-    # # *************************************************************************
-
-    # # there can only one outgoing arc at most, if there are no incoming arcs
-
-    # def rule_constr_max_outgoing_arcs(m,g,l):
-
-    #     # the number of predefined incoming arcs
-
-    #     n_in_pre = sum(
-    #         len(m.set_J_pre[(g,l_star,l)]) # = n_in_pre
-    #         for l_star in m.set_L[g]
-    #         if l_star not in m.set_L_exp[g]
-    #         if l_star != l
-    #         )
-
-    #     # if there is at least one predefined incoming arc, skip constraint
-
-    #     if n_in_pre >= 1:
-
-    #         return pyo.Constraint.Skip
-
-    #     # the number of non-predefined incoming arcs
-
-    #     n_in_opt = sum(
-    #         len(m.set_J_new[(g,l_star,l)]) # = n_in_pre
-    #         for l_star in m.set_L[g]
-    #         if l_star not in m.set_L_exp[g]
-    #         if l_star != l
-    #         )
-
-    #     n_in_max = n_in_pre + n_in_opt
-
-    #     # the number of predefined outgoing arcs
-
-    #     n_out_pre = sum(
-    #         len(m.set_J_pre[(g,l,l_line)])
-    #         for l_line in m.set_L[g]
-    #         if l_line not in m.set_L_imp[g]
-    #         if l_line != l
-    #         )
-
-    #     # the constraint is infeasible if the maximum number of incoming arcs
-    #     # is zero and the number of predefined outgoing arcs is bigger than 1
-
-    #     if n_in_max == 0 and n_out_pre >= 2:
-
-    #         return pyo.Constraint.Infeasible
-
-    #     # DONE: it is also infeasible if the maximum number of incoming arcs is
-    #     # zero and the number of predefined outgoing arcs is one and the poten-
-    #     # tial outgoing arcs include mandatory arcs (i.e. sum(...)=1 )
-
-    #     n_out_fcd = sum(
-    #         len(m.set_J_mdt[(g,l,l_line)])
-    #         for l_line in m.set_L[g]
-    #         if l_line not in m.set_L_imp[g]
-    #         if l_line != l
-    #         )
-
-    #     if n_in_max == 0 and n_out_pre == 1 and n_out_fcd >= 1:
-
-    #         return pyo.Constraint.Infeasible
-
-    #     # the number of non-predefined outgoing arcs
-
-    #     n_out_opt = sum(
-    #         len(m.set_J_new[(g,l,l_line)])
-    #         for l_line in m.set_L[g]
-    #         if l_line not in m.set_L_imp[g]
-    #         if l_line != l
-    #         )
-
-    #     n_out_max = n_out_pre + n_out_opt
-
-    #     if n_out_max <= 1:
-
-    #         # there can only be one outgoing arc at most: redundant constraint
-
-    #         return pyo.Constraint.Skip
-
-    #     else: # more than one outgoing arc is possible
-
-    #         M_gl = n_out_max - 1
-
-    #         return (
-    #             sum(
-    #                 sum(
-    #                     sum(m.var_delta_arc_inv_glljh[(g,l,l_diamond,j,s)]
-    #                         for s in m.set_H_gllj[(g,l,l_diamond,j)])
-    #                     for j in m.set_J_new[(g,l,l_diamond)]
-    #                     )
-    #                 #+len(m.set_J_pre[(g,l,l_diamond)]) # = n_out_pre
-    #                 for l_diamond in m.set_L[g]
-    #                 if l_diamond not in m.set_L_imp[g]
-    #                 if l_diamond != l
-    #                 )+n_out_pre
-    #             <= 1 + M_gl*
-    #             sum(
-    #                 sum(
-    #                     sum(m.var_delta_arc_inv_glljh[
-    #                         (g,l_star,l,j_star,s_star)]
-    #                         for s_star in m.set_H_gllj[(g,l_star,l,j_star)])
-    #                     for j_star in m.set_J_new[(g,l_star,l)]
-    #                     )
-    #                 #+len(m.set_J_pre[(g,l_star,l)]) # = n_in_pre
-    #                 for l_star in m.set_L[g]
-    #                 if l_star not in m.set_L_exp[g]
-    #                 if l_star != l
-    #                 )+n_in_pre
-    #             )
-
-    # model.constr_max_outgoing_arcs = pyo.Constraint(
-    #     model.set_GL_not_exp_imp,
-    #     rule=rule_constr_max_outgoing_arcs)
+    
+    add_network_restrictions(model)
 
     # *************************************************************************
     # *************************************************************************
diff --git a/tests/test_esipp_network.py b/tests/test_esipp_network.py
index 2617dbf..d731bb9 100644
--- a/tests/test_esipp_network.py
+++ b/tests/test_esipp_network.py
@@ -2170,12 +2170,10 @@ class TestNetwork:
     # *************************************************************************
 
     def test_tree_topology(self):
+        
         # create a network object with a tree topology
-
         tree_network = binomial_tree(3, create_using=MultiDiGraph)
-
-        network = Network(tree_network)
-
+        network = Network(incoming_graph_data=tree_network)
         for edge_key in network.edges(keys=True):
             arc = ArcsWithoutLosses(
                 name=str(edge_key),
@@ -2184,44 +2182,36 @@ class TestNetwork:
                 specific_capacity_cost=0,
                 capacity_is_instantaneous=False,
             )
-
             network.add_edge(*edge_key, **{Network.KEY_ARC_TECH: arc})
-
+            
         # assert that it does not have a tree topology
-
         assert not network.has_tree_topology()
 
         # select all the nodes
-
         for edge_key in network.edges(keys=True):
             network.edges[edge_key][Network.KEY_ARC_TECH].options_selected[0] = True
-
+            
         # assert that it has a tree topology
-
         assert network.has_tree_topology()
 
     # *************************************************************************
     # *************************************************************************
 
     def test_pseudo_unique_key_generation(self):
+        
         # create network
-
         network = Network()
 
         # add node A
-
         network.add_waypoint_node(node_key="A")
 
         # add node B
-
         network.add_waypoint_node(node_key="B")
 
         # identify nodes
-
         network.identify_node_types()
 
         # add arcs
-
         key_list = [
             "3e225573-4e78-48c8-bb08-efbeeb795c22",
             "f6d30428-15d1-41e9-a952-0742eaaa5a31",
diff --git a/tests/test_esipp_problem.py b/tests/test_esipp_problem.py
index 8e1b32a..8cf5126 100644
--- a/tests/test_esipp_problem.py
+++ b/tests/test_esipp_problem.py
@@ -226,7 +226,7 @@ class TestESIPPProblem:
         # ipp.instantiate(place_fixed_losses_upstream_if_possible=False)
 
         ipp.instantiate(initialise_ancillary_sets=init_aux_sets)
-        ipp.instance.pprint()
+        
         # optimise
         ipp.optimise(
             solver_name=solver,
@@ -5996,7 +5996,7 @@ class TestESIPPProblem:
                 solver_options={},
                 perform_analysis=False,
                 plot_results=False,  
-                print_solver_output=True,
+                print_solver_output=False,
                 time_frame=tf,
                 networks={"mynet": mynet},
                 static_losses_mode=static_losses_mode,
@@ -6147,7 +6147,7 @@ class TestESIPPProblem:
                 solver_options={},
                 perform_analysis=False,
                 plot_results=False,  # True,
-                print_solver_output=True,
+                print_solver_output=False,
                 time_frame=tf,
                 networks={"mynet": mynet},
                 static_losses_mode=static_losses_mode,
-- 
GitLab