How to use the pyomo.environ.Set function in Pyomo

def build(self):
        super(_PhysicalParameterBlock, self).build()

        self.p1 = Phase()
        self.p2 = Phase()

        self.c1 = Component()
        self.c2 = Component()

        self.phase_equilibrium_idx = Set(initialize=["e1", "e2"])
        self.element_list = Set(initialize=["H", "He", "Li"])
        self.element_comp = {"c1": {"H": 1, "He": 2, "Li": 3},
                             "c2": {"H": 4, "He": 5, "Li": 6}}

        self.phase_equilibrium_list = \
            {"e1": ["c1", ("p1", "p2")],
             "e2": ["c2", ("p1", "p2")]}

        # Attribute to switch flow basis for testing
        self.basis_switch = 1
        self.default_balance_switch = 1

        self._state_block_class = TestStateBlock

for (i, j) in model.Arcs:
            if i == node:
                retval.append(j)
        return retval

    model.NodesOut = py.Set(model.Nodes, initialize=nodes_out_init)

    # compute NodesIn, i.e., set of nodes connected to considered node by ingoing arc
    def nodes_in_init(model, node):
        retval = []
        for (i, j) in model.Arcs:
            if j == node:
                retval.append(i)
        return retval

    model.NodesIn = py.Set(model.Nodes, initialize=nodes_in_init)

    # add flow balance constraints corresponding to the node demands
    def flow_balance_rule(model, node):
        return sum(model.Flow[i, node] for i in model.NodesIn[node]) \
               - sum(model.Flow[node, j] for j in model.NodesOut[node]) \
               == model.Demand[node]

    model.FlowBalance_cons = py.Constraint(model.Nodes, rule=flow_balance_rule)

    # compute unique flow-path P(startNode,endNode) from entry to exit; given by list of nodes of the path
    pathNodes = nx.shortest_path(graph, source=startNode, target=endNode)
    # non zero coefficients of objective function
    dic_arc_coef = {}
    # determine coefficients for objective function
    # if for an arc (u,v), u, respectively v, are not in pathNodes, then the coefficient is 0
    # if arc (u,v) of pathNodes satisfies P(startNode, u) subset P(startNode,v), then coefficient is 1, otherwise -1

return dic_nodes_MinCapacity[node], dic_nodes_MaxCapacity[node]

        model.Demand = py.Var(model.Nodes, bounds=demandCapacities)

    # create arc flow variables for every arc of the network
    model.Flow = py.Var(model.Arcs)

    # compute NodesOut, i.e., set of nodes that are connected to considered node by outgoing arc
    def nodes_out_init(model, node):
        retval = []
        for (i, j) in model.Arcs:
            if i == node:
                retval.append(j)
        return retval

    model.NodesOut = py.Set(model.Nodes, initialize=nodes_out_init)

    # compute NodesIn, i.e., set of nodes connected to considered node by ingoing arc
    def nodes_in_init(model, node):
        retval = []
        for (i, j) in model.Arcs:
            if j == node:
                retval.append(i)
        return retval

    model.NodesIn = py.Set(model.Nodes, initialize=nodes_in_init)

    # add flow balance constraints corresponding to the node demands
    def flow_balance_rule(model, node):
        return sum(model.Flow[i, node] for i in model.NodesIn[node]) \
               - sum(model.Flow[node, j] for j in model.NodesOut[node]) \
               == model.Demand[node]

def build(self):
        '''
        Callable method for Block construction.
        '''
        super(BTParameterData, self).build()

        self.cubic_type = CubicEoS.PR

        # Add Component objects
        self.benzene = Component()
        self.toluene = Component()

        # List of phase equilibrium index
        self.phase_equilibrium_idx = Set(initialize=[1, 2])

        self.phase_equilibrium_list = \
            {1: ["benzene", ("Vap", "Liq")],
             2: ["toluene", ("Vap", "Liq")]}

        # Thermodynamic reference state
        self.pressure_ref = Param(mutable=True,
                                  default=101325,
                                  doc='Reference pressure [Pa]')
        self.temperature_ref = Param(mutable=True,
                                     default=298.15,
                                     doc='Reference temperature [K]')

        # Critical Properties
        pressure_crit_data = {'benzene': 48.9e5,
                              'toluene': 41.0e5}

def rule_pmax(m, n):
        if n == 'n13':
            return 41.0
        return 70.0
    model.pmax = aml.Param(model.NODE, initialize=rule_pmax, mutable=True)

    # supply
    model.SUP = aml.Set(initialize=[1])
    model.sloc = aml.Param(model.SUP, initialize='n1')
    model.smin = aml.Param(model.SUP, within=aml.NonNegativeReals, initialize=0.000, mutable=True)
    model.smax = aml.Param(model.SUP, within=aml.NonNegativeReals, initialize=30, mutable=True)
    model.scost = aml.Param(model.SUP, within=aml.NonNegativeReals)

    # demand
    model.DEM = aml.Set(initialize=[1])
    model.dloc = aml.Param(model.DEM, initialize='n13')
    model.d = aml.Param(model.DEM, within=aml.PositiveReals, initialize=10, mutable=True)

    # physical data

    model.TDEC = aml.Param(initialize=9.5)

    model.eps = aml.Param(initialize=0.025, within=aml.PositiveReals)
    model.z = aml.Param(initialize=0.80, within=aml.PositiveReals)
    model.rhon = aml.Param(initialize=0.72, within=aml.PositiveReals)
    model.R = aml.Param(initialize=8314.0, within=aml.PositiveReals)
    model.M = aml.Param(initialize=18.0, within=aml.PositiveReals)
    model.pi = aml.Param(initialize=3.14, within=aml.PositiveReals)
    model.nu2 = aml.Param(within=aml.PositiveReals,mutable=True)
    model.lam = aml.Param(model.LINK, within=aml.PositiveReals, mutable=True)
    model.A = aml.Param(model.LINK, within=aml.NonNegativeReals, mutable=True)

def remove_relaxation(self):
        """
        Remove any auto-created vars/constraints from the relaxation block
        """
        # this default implementation should work for most relaxations
        # it removes all vars and constraints on this block data object
        self._remove_from_persistent_solvers()
        comps = [pe.Block, pe.Constraint, pe.Var, pe.Set, pe.Param]
        for comp in comps:
            comps_to_del = list(self.component_objects([comp], descend_into=False))
            for _comp in comps_to_del:
                self.del_component(_comp)
        for comp in comps:
            comps_to_del = list(self.component_data_objects([comp], descend_into=False))
            for _comp in comps_to_del:
                self.del_component(_comp)

def build(self):
        '''
        Callable method for Block construction.
        '''
        super(HDAReactionParameterData, self).build()

        self.reaction_block_class = HDAReactionBlock

        # List of valid phases in property package
        self.phase_list = Set(initialize=['Vap'])

        # Component list - a list of component identifiers
        self.component_list = Set(initialize=['benzene',
                                              'toluene',
                                              'hydrogen',
                                              'methane'])

        # Reaction Index
        self.rate_reaction_idx = Set(initialize=["R1"])

        # Reaction Stoichiometry
        self.rate_reaction_stoichiometry = {("R1", "Vap", "benzene"): 1,
                                            ("R1", "Vap", "toluene"): -1,
                                            ("R1", "Vap", "hydrogen"): -1,
                                            ("R1", "Vap", "methane"): 1,
                                            ("R1", "Liq", "benzene"): 0,

yub = pyo.value(y.ub)

    check_var_pts(x, x_pts=x_pts)
    check_var_pts(y)

    if x.is_fixed() and y.is_fixed():
        b.xy_fixed_eq = pyo.Constraint(expr= w == pyo.value(x) * pyo.value(y))
    elif x.is_fixed():
        b.x_fixed_eq = pyo.Constraint(expr= w == pyo.value(x) * y)
    elif y.is_fixed():
        b.y_fixed_eq = pyo.Constraint(expr= w == x * pyo.value(y))
    elif len(x_pts) == 2:
        _build_mccormick_relaxation(b, x=x, y=y, w=w, relaxation_side=relaxation_side)
    else:
        # create the lambda variables (binaries for the pw representation)
        b.interval_set = pyo.Set(initialize=range(1, len(x_pts)))
        b.lam = pyo.Var(b.interval_set, within=pyo.Binary)

        # create the delta y variables
        b.delta_y = pyo.Var(b.interval_set, bounds=(0, None))

        # create the "sos1" constraint
        b.lam_sos1 = pyo.Constraint(expr=sum(b.lam[n] for n in b.interval_set) == 1.0)

        # create the x interval constraints
        b.x_interval_lb = pyo.Constraint(expr=sum(x_pts[n - 1] * b.lam[n] for n in b.interval_set) <= x)
        b.x_interval_ub = pyo.Constraint(expr=x <= sum(x_pts[n] * b.lam[n] for n in b.interval_set))

        # create the y constraints
        b.y_con = pyo.Constraint(expr=y == ylb + sum(b.delta_y[n] for n in b.interval_set))

        def delta_yn_ub_rule(m, n):