src/topupopt/problems/esipp/blocks/converters.py

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

src/topupopt/problems/esipp/network.py

@@ -612,21 +612,62 @@ class Network(nx.MultiDiGraph):
         KEY_ARC_TECH_CAPACITY_INSTANTANEOUS,
         KEY_ARC_TECH_STATIC_LOSS,
     )
+    
+    NET_TYPE_HYBRID = 0
+    NET_TYPE_TREE = 1
+    NET_TYPE_REV_TREE = 2
+    
+    NET_TYPES = (
+        NET_TYPE_HYBRID,
+        NET_TYPE_TREE,
+        NET_TYPE_REV_TREE
+        )
 
-    def __init__(self, incoming_graph_data=None, **attr):
+    def __init__(self, network_type = NET_TYPE_HYBRID, **kwargs):
         # run base class init routine
 
-        nx.MultiDiGraph.__init__(self, incoming_graph_data=incoming_graph_data, **attr)
+        nx.MultiDiGraph.__init__(self, **kwargs)
 
         # identify node types
 
         self.identify_node_types()
 
         # declare variables for the nodes without directed arc limitations
+        self.network_type = network_type
+
+        self.nodes_w_in_dir_arc_limitations = dict()
+    
+        self.nodes_w_out_dir_arc_limitations = dict()
+
+    # *************************************************************************
+    # *************************************************************************
+    
+    def _set_up_node(self, node_key, max_number_in_arcs: int = None, max_number_out_arcs: int = None):
+        
+        if self.should_be_tree_network():
+            # nodes have to be part of a tree: one incoming arc per node at most
+            self.nodes_w_in_dir_arc_limitations[node_key] = 1
+        elif self.should_be_reverse_tree_network():
+            # nodes have to be part of a reverse tree: one outgoing arc per node at most
+            self.nodes_w_out_dir_arc_limitations[node_key] = 1
+        else: 
+            # nodes have no peculiar restrictions or they are defined 1 by 1
+            if type(max_number_in_arcs) != type(None):
+                self.nodes_w_in_dir_arc_limitations[node_key] = max_number_in_arcs
+            if type(max_number_out_arcs) != type(None):
+                self.nodes_w_out_dir_arc_limitations[node_key] = max_number_out_arcs
 
-        self.nodes_wo_in_dir_arc_limitations = []
+    # *************************************************************************
+    # *************************************************************************
+    
+    def should_be_tree_network(self) -> bool:
+        return self.network_type == self.NET_TYPE_TREE
 
-        self.nodes_wo_out_dir_arc_limitations = []
+    # *************************************************************************
+    # *************************************************************************
+    
+    def should_be_reverse_tree_network(self) -> bool:
+        return self.network_type == self.NET_TYPE_REV_TREE
 
     # *************************************************************************
     # *************************************************************************
@@ -661,23 +702,25 @@ class Network(nx.MultiDiGraph):
 
     # add a new supply/demand node
 
-    def add_source_sink_node(self, node_key, base_flow: dict):
+    def add_source_sink_node(self, node_key, base_flow: dict, **kwargs):
         node_dict = {
             self.KEY_NODE_TYPE: self.KEY_NODE_TYPE_SOURCE_SINK,
             self.KEY_NODE_BASE_FLOW: base_flow,
         }
 
         self.add_node(node_key, **node_dict)
+        self._set_up_node(node_key, **kwargs)
 
     # *************************************************************************
     # *************************************************************************
 
     # add a new waypoint node
 
-    def add_waypoint_node(self, node_key):
+    def add_waypoint_node(self, node_key, **kwargs):
         node_dict = {self.KEY_NODE_TYPE: self.KEY_NODE_TYPE_WAY}
 
         self.add_node(node_key, **node_dict)
+        self._set_up_node(node_key, **kwargs)
 
     # *************************************************************************
     # *************************************************************************

src/topupopt/problems/esipp/problem.py

@@ -63,6 +63,15 @@ class InfrastructurePlanningProblem(EnergySystem):
         STATIC_LOSS_MODE_US,
         STATIC_LOSS_MODE_DS,
     )
+    
+    NODE_PRICE_LAMBDA = 1
+    NODE_PRICE_DELTA = 2
+    NODE_PRICE_OTHER = 3
+    NODE_PRICES = (
+        NODE_PRICE_LAMBDA,
+        NODE_PRICE_DELTA,
+        NODE_PRICE_OTHER
+        )
 
     # *************************************************************************
     # *************************************************************************
@@ -80,6 +89,7 @@ class InfrastructurePlanningProblem(EnergySystem):
         converters: dict = None,
         prepare_model: bool = True,
         validate_inputs: bool = True,
+        node_price_model = NODE_PRICE_DELTA
     ):  # TODO: switch to False when everything is more mature
         # *********************************************************************
 
@@ -1830,22 +1840,14 @@ class InfrastructurePlanningProblem(EnergySystem):
         }
 
         set_L_max_in_g = {
-            g: tuple(
-                l
-                for l in self.networks[g].nodes
-                if l not in self.networks[g].nodes_wo_in_dir_arc_limitations
-            )
+            g: tuple(self.networks[g].nodes_w_in_dir_arc_limitations.keys())
             for g in self.networks.keys()
-        }
+            }
 
-        # set_L_max_out_g = {
-        #     g: tuple(
-        #         l
-        #         for l in self.networks[g].nodes
-        #         if l not in self.networks[g].nodes_wo_out_dir_arc_limitations
-        #         )
-        #     for g in self.networks.keys()
-        #     }
+        set_L_max_out_g = {
+            g: tuple(self.networks[g].nodes_w_out_dir_arc_limitations.keys())
+            for g in self.networks.keys()
+            }
 
         set_GL = tuple((g, l) for g in set_G for l in set_L[g])
 
@@ -1897,7 +1899,7 @@ class InfrastructurePlanningProblem(EnergySystem):
             for (g, l) in set_GL_exp_imp
             for (q, p, k) in set_QPK
         }
-
+        
         # set of GLKS tuples
         set_GLQPKS = tuple(
             (*glqpk, s) for glqpk, s_tuple in set_S.items() for s in s_tuple
@@ -2547,6 +2549,17 @@ class InfrastructurePlanningProblem(EnergySystem):
             for s in set_S[(g, l, q, p, k)]
         }
 
+        # price function convexity
+
+        param_price_function_is_convex = {
+            (g, l, q, p, k): (
+                self.networks[g].nodes[l][Network.KEY_NODE_PRICES][(q, p, k)].price_monotonically_increasing_with_volume()
+                if l in set_L_imp[g] else
+                self.networks[g].nodes[l][Network.KEY_NODE_PRICES][(q, p, k)].price_monotonically_decreasing_with_volume()
+                )
+            for (g, l, q, p, k) in set_S
+        }
+
         # maximum resource volume per segment (infinity is the default)
 
         param_v_max_glqpks = {
@@ -3317,7 +3330,7 @@ class InfrastructurePlanningProblem(EnergySystem):
                 "set_L_imp": set_L_imp,
                 "set_L_exp": set_L_exp,
                 "set_L_max_in_g": set_L_max_in_g,
-                #'set_L_max_out_g': set_L_max_out_g,
+                'set_L_max_out_g': set_L_max_out_g,
                 "set_GL": set_GL,
                 "set_GL_exp": set_GL_exp,
                 "set_GL_imp": set_GL_imp,
@@ -3449,6 +3462,7 @@ class InfrastructurePlanningProblem(EnergySystem):
                 "param_c_df_qp": param_c_df_qp,
                 "param_c_time_qpk": param_c_time_qpk,
                 "param_p_glqpks": param_p_glqpks,
+                "param_price_function_is_convex": param_price_function_is_convex,
                 "param_v_max_glqpks": param_v_max_glqpks,
                 # *****************************************************************
                 # converters

src/topupopt/problems/esipp/resource.py

@@ -12,7 +12,11 @@ from numbers import Real
 class ResourcePrice:
     """A class for piece-wise linear resource prices in network problems."""
 
-    def __init__(self, prices: list or int, volumes: list = None):
+    def __init__(
+            self, 
+            prices: list or int, 
+            volumes: list = None
+            ):
         # how do we keep the size of the object as small as possible
         # if the tariff is time-invariant, how can information be stored?
         # - a flag
@@ -206,30 +210,10 @@ class ResourcePrice:
 
     # *************************************************************************
     # *************************************************************************
-
-    def is_equivalent(self, other) -> bool:
-        """Returns True if a given ResourcePrice is equivalent to another."""
-        # resources are equivalent if:
-        # 1) the prices are the same
-        # 2) the volume limits are the same
-
-        # the number of segments has to match
-        if self.number_segments() != other.number_segments():
-            return False  # the number of segments do not match
-        # check the prices
-        if self.prices != other.prices:
-            return False  # prices are different
-        # prices match, check the volumes
-        if self.volumes != other.volumes:
-            return False  # volumes are different
-        return True  # all conditions have been met
-
-    # *************************************************************************
-    # *************************************************************************
     
     def __eq__(self, o) -> bool:
         """Returns True if a given ResourcePrice is equivalent to another."""
-        return self.is_equivalent(o)
+        return hash(self) == hash(o)
     
     def __hash__(self):
         return hash(
@@ -260,9 +244,7 @@ def are_prices_time_invariant(resource_prices_qpk: dict) -> bool:
     # check if the tariffs per period and assessment are equivalent
     for qp, qpk_list in qpk_qp.items():
         for i in range(len(qpk_list) - 1):
-            if not resource_prices_qpk[qpk_list[0]].is_equivalent(
-                resource_prices_qpk[qpk_list[i + 1]]
-            ):
+            if not resource_prices_qpk[qpk_list[0]] == resource_prices_qpk[qpk_list[i + 1]]:
                 return False
     # all tariffs are equivalent per period and assessment: they are invariant
     return True

src/topupopt/problems/esipp/utils.py

@@ -99,7 +99,7 @@ def statistics(ipp: InfrastructurePlanningProblem,
     imports_qpk = {
         qpk: pyo.value( 
             sum(
-                ipp.instance.var_if_glqpks[(g,l_imp,*qpk, s)]
+                ipp.instance.var_trans_flows_glqpks[(g,l_imp,*qpk, s)]
                 for g, l_imp in import_node_keys
                 # for g in ipp.networks
                 # for l_imp in ipp.networks[g].import_nodes
@@ -114,7 +114,7 @@ def statistics(ipp: InfrastructurePlanningProblem,
     exports_qpk = {
         qpk: pyo.value(
             sum(
-                ipp.instance.var_ef_glqpks[(g,l_exp,*qpk, s)]
+                ipp.instance.var_trans_flows_glqpks[(g,l_exp,*qpk, s)]
                 for g, l_exp in export_node_keys
                 # for g in ipp.networks
                 # for l_exp in ipp.networks[g].export_nodes

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

tests/test_esipp_resource.py

@@ -132,8 +132,8 @@ class TestResourcePrice:
         volumes = None
         res_p1 = ResourcePrice(prices=prices, volumes=volumes)
         res_p2 = ResourcePrice(prices=[prices], volumes=[volumes])
-        assert res_p1.is_equivalent(res_p2)
-        assert res_p2.is_equivalent(res_p1)
+        assert res_p1 == res_p2
+        assert res_p2 == res_p1
 
         # *********************************************************************
 
@@ -144,8 +144,8 @@ class TestResourcePrice:
         volumes = None
         res_p1 = ResourcePrice(prices=prices, volumes=volumes)
         res_p2 = ResourcePrice(prices=[prices + 1], volumes=[volumes])
-        assert not res_p1.is_equivalent(res_p2)
-        assert not res_p2.is_equivalent(res_p1)
+        assert not res_p1 == res_p2
+        assert not res_p2 == res_p1
 
         # *********************************************************************
 
@@ -156,8 +156,8 @@ class TestResourcePrice:
         volumes = None
         res_p1 = ResourcePrice(prices=prices, volumes=volumes)
         res_p2 = ResourcePrice(prices=prices, volumes=volumes)
-        assert res_p1.is_equivalent(res_p2)
-        assert res_p2.is_equivalent(res_p1)
+        assert res_p1 == res_p2
+        assert res_p2 == res_p1
 
         # *********************************************************************
 
@@ -168,8 +168,8 @@ class TestResourcePrice:
         volumes = None
         res_p1 = ResourcePrice(prices=prices, volumes=volumes)
         res_p2 = ResourcePrice(prices=prices + 1, volumes=volumes)
-        assert not res_p1.is_equivalent(res_p2)
-        assert not res_p2.is_equivalent(res_p1)
+        assert not res_p1 == res_p2
+        assert not res_p2 == res_p1
 
         # *********************************************************************
         # *********************************************************************
@@ -183,8 +183,8 @@ class TestResourcePrice:
         volumes = 1
         res_p1 = ResourcePrice(prices=prices, volumes=volumes)
         res_p2 = ResourcePrice(prices=[prices], volumes=[volumes])
-        assert res_p1.is_equivalent(res_p2)
-        assert res_p2.is_equivalent(res_p1)
+        assert res_p1 == res_p2
+        assert res_p2 == res_p1
 
         # *********************************************************************
 
@@ -195,8 +195,8 @@ class TestResourcePrice:
         volumes = 1
         res_p1 = ResourcePrice(prices=prices, volumes=volumes)
         res_p2 = ResourcePrice(prices=[prices + 1], volumes=[volumes])
-        assert not res_p1.is_equivalent(res_p2)
-        assert not res_p2.is_equivalent(res_p1)
+        assert not res_p1 == res_p2
+        assert not res_p2 == res_p1
 
         # *********************************************************************
 
@@ -207,8 +207,8 @@ class TestResourcePrice:
         volumes = 1
         res_p1 = ResourcePrice(prices=prices, volumes=volumes)
         res_p2 = ResourcePrice(prices=prices, volumes=volumes)
-        assert res_p1.is_equivalent(res_p2)
-        assert res_p2.is_equivalent(res_p1)
+        assert res_p1 == res_p2
+        assert res_p2 == res_p1
 
         # *********************************************************************
 
@@ -219,8 +219,8 @@ class TestResourcePrice:
         volumes = 1
         res_p1 = ResourcePrice(prices=prices, volumes=volumes)
         res_p2 = ResourcePrice(prices=prices + 1, volumes=volumes)
-        assert not res_p1.is_equivalent(res_p2)
-        assert not res_p2.is_equivalent(res_p1)
+        assert not res_p1 == res_p2
+        assert not res_p2 == res_p1
 
         # *********************************************************************
 
@@ -231,8 +231,8 @@ class TestResourcePrice:
         volumes = 1
         res_p1 = ResourcePrice(prices=prices, volumes=volumes)
         res_p2 = ResourcePrice(prices=[prices], volumes=[volumes + 1])
-        assert not res_p1.is_equivalent(res_p2)
-        assert not res_p2.is_equivalent(res_p1)
+        assert not res_p1 == res_p2
+        assert not res_p2 == res_p1
 
         # *********************************************************************
 
@@ -243,8 +243,8 @@ class TestResourcePrice:
         volumes = 1
         res_p1 = ResourcePrice(prices=prices, volumes=volumes)
         res_p2 = ResourcePrice(prices=prices, volumes=volumes + 1)
-        assert not res_p1.is_equivalent(res_p2)
-        assert not res_p2.is_equivalent(res_p1)
+        assert not res_p1 == res_p2
+        assert not res_p2 == res_p1
 
         # *********************************************************************
 
@@ -255,8 +255,8 @@ class TestResourcePrice:
         volumes = 1
         res_p1 = ResourcePrice(prices=prices, volumes=volumes)
         res_p2 = ResourcePrice(prices=[prices], volumes=[None])
-        assert not res_p1.is_equivalent(res_p2)
-        assert not res_p2.is_equivalent(res_p1)
+        assert not res_p1 == res_p2
+        assert not res_p2 == res_p1
 
         # *********************************************************************
 
@@ -267,8 +267,8 @@ class TestResourcePrice:
         volumes = 1
         res_p1 = ResourcePrice(prices=prices, volumes=volumes)
         res_p2 = ResourcePrice(prices=prices, volumes=None)
-        assert not res_p1.is_equivalent(res_p2)
-        assert not res_p2.is_equivalent(res_p1)
+        assert not res_p1 == res_p2
+        assert not res_p2 == res_p1
 
         # *********************************************************************
         # *********************************************************************
@@ -294,8 +294,8 @@ class TestResourcePrice:
         volumes = [1, None]
         res_p1 = ResourcePrice(prices=prices, volumes=volumes)
         res_p2 = ResourcePrice(prices=prices, volumes=volumes)
-        assert res_p1.is_equivalent(res_p2)
-        assert res_p2.is_equivalent(res_p1)
+        assert res_p1 == res_p2
+        assert res_p2 == res_p1
 
         # two segments, no volume limit, same format
         # prices do not match = False
@@ -306,8 +306,8 @@ class TestResourcePrice:
         prices = [2, 3]
         volumes = [1, None]
         res_p2 = ResourcePrice(prices=prices, volumes=volumes)
-        assert not res_p1.is_equivalent(res_p2)
-        assert not res_p2.is_equivalent(res_p1)
+        assert not res_p1 == res_p2
+        assert not res_p2 == res_p1
 
         # *********************************************************************
 
@@ -320,8 +320,8 @@ class TestResourcePrice:
         volumes = [1, 3]
         res_p1 = ResourcePrice(prices=prices, volumes=volumes)
         res_p2 = ResourcePrice(prices=prices, volumes=volumes)
-        assert res_p1.is_equivalent(res_p2)
-        assert res_p2.is_equivalent(res_p1)
+        assert res_p1 == res_p2
+        assert res_p2 == res_p1
 
         # two segments, volume limit, same format: False
         # prices do not match = False
@@ -332,8 +332,8 @@ class TestResourcePrice:
         prices = [1, 4]
         volumes = [1, 4]
         res_p2 = ResourcePrice(prices=prices, volumes=volumes)
-        assert not res_p1.is_equivalent(res_p2)
-        assert not res_p2.is_equivalent(res_p1)
+        assert not res_p1 == res_p2
+        assert not res_p2 == res_p1
 
         # *********************************************************************
 
@@ -348,8 +348,8 @@ class TestResourcePrice:
         prices = [1, 3]
         volumes = [1, 5]
         res_p2 = ResourcePrice(prices=prices, volumes=volumes)
-        assert not res_p1.is_equivalent(res_p2)
-        assert not res_p2.is_equivalent(res_p1)
+        assert not res_p1 == res_p2
+        assert not res_p2 == res_p1
 
         # single segment, volume limit, same format
         # volumes do not match = False
@@ -360,8 +360,8 @@ class TestResourcePrice:
         prices = [1, 3]
         volumes = [1, None]
         res_p2 = ResourcePrice(prices=prices, volumes=volumes)
-        assert not res_p1.is_equivalent(res_p2)
-        assert not res_p2.is_equivalent(res_p1)
+        assert not res_p1 == res_p2
+        assert not res_p2 == res_p1
 
         # *********************************************************************
         # *********************************************************************
@@ -374,8 +374,8 @@ class TestResourcePrice:
         prices = [1, 3, 5]
         volumes = [1, 4, None]
         res_p2 = ResourcePrice(prices=prices, volumes=volumes)
-        assert not res_p1.is_equivalent(res_p2)
-        assert not res_p2.is_equivalent(res_p1)
+        assert not res_p1 == res_p2
+        assert not res_p2 == res_p1
 
         # *********************************************************************
         # *********************************************************************