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

def fr1(cracmodel):
            # A <> B + H
            return cracmodel.r1 == k[0] * pyo.exp(-(E[0]/(gc))*((1/T)-(1/Tr))) * (cracmodel.a - (cracmodel.b * cracmodel.h)/keq(T)) * (1.0/((1+kst*cracmodel.b)*(1+kco*cracmodel.i)))
        def fr2(cracmodel):
            # A > C + D
            return cracmodel.r2 == k[1] * pyo.exp(-(E[1]/(gc))*((1/T)-(1/Tr)))  * cracmodel.a * (1.0/(1+kco*cracmodel.i))

        def fr4(cracmodel):
            # D + 2H > 2F
            return cracmodel.r4 == k[2] * pyo.exp(-(E[2]/(gc))*((1/T)-(1/Tr)))  * cracmodel.d * cracmodel.h
        def fr5(cracmodel):
            # F+G > I + 4H
            return cracmodel.r5 == k[3] * pyo.exp(-(E[3]/(gc))*((1/T)-(1/Tr)))  * cracmodel.f * cracmodel.g

        cracmodel.er1 = pyo.Constraint( rule = fr1)
        cracmodel.er2 = pyo.Constraint( rule = fr2)
        cracmodel.er4 = pyo.Constraint( rule = fr4)
        cracmodel.er5 = pyo.Constraint( rule = fr5)
        cracmodel.sets = pyo.RangeSet(9)
        cracmodel.dum = pyo.Var(cracmodel.sets, domain = pyo.Reals)

        def fra(cracmodel):  # A - 1,2,3
            return cracmodel.dum[1]  == ca0-cracmodel.a - cracmodel.r1 -cracmodel.r2 
        def frb(cracmodel): # B - 1
            return cracmodel.dum[2]  == cb0-cracmodel.b+cracmodel.r1
        def frc(cracmodel): # C - 2
            return cracmodel.dum[3]  == cc0-cracmodel.c+cracmodel.r2
        def frd(cracmodel): # D - 2,4
            return cracmodel.dum[4]  == cd0-cracmodel.d+cracmodel.r2-cracmodel.r4
        def frf(cracmodel): # F - 3,4,5
            return cracmodel.dum[5]  == cf0-cracmodel.f+2*cracmodel.r4-cracmodel.r5

def rule_teq(b):
            return b._teq == 0.5 * (b._t1 + b.temperature_dew -
                                    sqrt((b._t1 - b.temperature_dew)**2 +
                                         b.eps_2**2))
        self._teq_constraint = Constraint(rule=rule_teq)

        def rule_tr_eq(b, i):
            return b._teq / b._params.temperature_crit[i]
        self._tr_eq = Expression(self._params.component_list,
                                 rule=rule_tr_eq,
                                 doc='Component reduced temperatures [-]')

        def rule_equilibrium(b, i):
            return b.fug_vap[i] == b.fug_liq[i]
        self.equilibrium_constraint = \
            Constraint(self._params.component_list, rule=rule_equilibrium)

def _enth_mol_phase_comp(self):
        self.enth_mol_phase_comp = Var(self._params.phase_list,
                                       self._params.component_list,
                                       doc="Phase-component molar specific "
                                           "enthalpies [J/mol]")

        def rule_enth_mol_phase_comp(b, p, j):
            if p == 'Vap':
                return b._enth_mol_comp_vap(j)
            else:
                return b._enth_mol_comp_liq(j)
        self.eq_enth_mol_phase_comp = Constraint(
            self._params.phase_list,
            self._params.component_list,
            rule=rule_enth_mol_phase_comp)

def declare_eq_branch_power_off(model, index_set, branches):
    """
    RHS of disjunct for declare_eq_branch_power_btheta_approx
    """
    m = model

    con_set = decl.declare_set("_con_eq_branch_power_off", model, index_set)

    m.eq_pf_branch_off = pe.Constraint(con_set)
    for branch_name in con_set:
        m.eq_pf_branch_off[branch_name] = \
            m.pf[branch_name] == 0.

# add inequality bound on objective functions if required
    # obj.expr <= objective_ub
    if objective_bound is not None:
        if len(deactivated_objectives) != 1:
            e = 'BoundsTightener: When providing objective_ub,' + \
                ' the model must have one and only one objective function.'
            logger.error(e)
            raise ValueError(e)
        original_obj = deactivated_objectives[0]
        if original_obj.sense == minimize:
            model.__objective_ineq = \
                pyo.Constraint(expr=original_obj.expr <= objective_bound)
        else:
            assert original_obj.sense == maximize
            model.__objective_ineq = pyo.Constraint(expr=original_obj.expr >= objective_bound)
        if isinstance(solver, PersistentSolver):
            solver.add_constraint(model.__objective_ineq)

    return initial_var_values, deactivated_objectives

initialize=101325,
                                  doc="vapor pressure ")

        self.x = Var(self.component_list_ref, initialize=1,
                     doc="temporary variable to compute vapor pressure")

        def rule_total_mass_balance(self):
            return self.flow_mol_phase['Liq'] + \
                self.flow_mol_phase['Vap'] == self.flow_mol
        self.eq_total = Constraint(rule=rule_total_mass_balance)

        def rule_comp_mass_balance(self, i):
            return self.flow_mol * self.mole_frac[i] == \
                self.flow_mol_phase['Liq'] * self.mole_frac_phase['Liq', i] + \
                self.flow_mol_phase['Vap'] * self.mole_frac_phase['Vap', i]
        self.eq_comp = Constraint(self.component_list_ref,
                                  rule=rule_comp_mass_balance)

        def rule_mole_frac(self):
            return sum(self.mole_frac_phase['Liq', i]
                       for i in self.component_list_ref) -\
                sum(self.mole_frac_phase['Vap', i]
                    for i in self.component_list_ref) == 0
        self.eq_sum_mol_frac = Constraint(rule=rule_mole_frac)

        if self.config.defined_state is False:
            # applied at outlet only
            self.eq_mol_frac_out = Constraint(expr=sum(self.mole_frac[i]
                                              for i in self.component_list_ref)
                                              == 1)

        if self.config.has_phase_equilibrium:

def _rate_constant(self):
        self.k_rxn = Var(initialize=1,
                         doc="Rate constant [m^3/mol.s]")

        try:
            self.arrhenius_eqn = Constraint(expr=self.k_rxn ==
                                            self.arrhenius_ref *
                                            exp(-self.energy_activation_ref /
                                                (self.gas_const_ref *
                                                 self.temperature_ref)))
        except AttributeError:
            # If constraint fails, clean up so that DAE can try again later
            self.del_component(self.k_rxn)
            self.del_component(self.arrhenius_eqn)
            raise

# PSE paper Eqn 13
            def rule_t1(b):
                return b._t1 == 0.5 * \
                    (b.temperature + b.temperature_bubble +
                     sqrt((b.temperature - b.temperature_bubble)**2 +
                          b.eps_1**2))
            self._t1_constraint = Constraint(rule=rule_t1)

            # PSE paper Eqn 14
            # TODO : Add option for supercritical extension
            def rule_teq(b):
                return b._teq == 0.5 * (b._t1 + b.temperature_dew -
                                        sqrt((b._t1 - b.temperature_dew)**2 +
                                             b.eps_2**2))
            self._teq_constraint = Constraint(rule=rule_teq)

            def rule_tr_eq(b, i):
                return b._teq / b._params.temperature_crit[i]
            self._tr_eq = Expression(self._params.component_list,
                                     rule=rule_tr_eq,
                                     doc='Component reduced temperatures [-]')

            def rule_equilibrium(b, i):
                return b.fug_vap[i] == b.fug_liq[i]
            self.equilibrium_constraint = \
                Constraint(self._params.component_list, rule=rule_equilibrium)