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

# heat capacity ratio
        self.heat_capacity_ratio = Expression(expr=self.cp_mol / self.cv_mol)
        self.scaling_factor[self.heat_capacity_ratio] = 1e1
        # Flows
        self.flow_vol = Expression(
            expr=self.flow_mol / self.dens_mol,
            doc="Total liquid + vapor volumetric flow (m3/s)",
        )
        self.scaling_factor[self.flow_vol] = 100

        self.flow_mass = Expression(
            expr=self.mw * self.flow_mol, doc="mass flow rate [kg/s]"
        )
        self.scaling_factor[self.flow_mass] = 1

        self.enth_mass = Expression(expr=self.enth_mol / mw, doc="Mass enthalpy (J/kg)")
        self.scaling_factor[self.enth_mass] = 1

        # Set the state vars dictionary
        if self.state_vars == StateVars.PH:
            self._state_vars_dict = {
                "flow_mol": self.flow_mol,
                "enth_mol": self.enth_mol,
                "pressure": self.pressure,
            }
        elif self.state_vars == StateVars.TPX and phase_set in (
            PhaseType.MIX,
            PhaseType.LG,
        ):
            self._state_vars_dict = {
                "flow_mol": self.flow_mol,
                "temperature": self.temperature,

model.asymptote = pyo.Var(initialize=15)
    model.rate_constant = pyo.Var(initialize=0.5)

    def response_rule(m):
        expr = m.asymptote * (1 - pyo.exp(-m.rate_constant * h))
        return expr

    model.response_function = pyo.Expression(rule=response_rule)

    def SSE_rule(m):
        return (y - m.response_function) ** 2

    def ComputeFirstStageCost_rule(m):
        return 0

    model.FirstStageCost = pyo.Expression(rule=ComputeFirstStageCost_rule)
    model.SecondStageCost = pyo.Expression(rule=SSE_rule)

    def total_cost_rule(model):
        return model.FirstStageCost + model.SecondStageCost

    # This objective is not needed by parmest
    model.Total_Cost_Objective = pyo.Objective(rule=total_cost_rule, 
                                               sense=pyo.minimize)
    
    return model

for j in b.params.component_list)
                for i in b.params.component_list)
        blk.am = Expression(blk.params.phase_list, rule=rule_am)

        def rule_am_eq(b, p):
            return sum(sum(
                b.mole_frac_phase_comp[p, i]*b.mole_frac_phase_comp[p, j] *
                sqrt(b._a_eq[i]*b._a_eq[j])*(1-b.params.kappa[i, j])
                for j in b.params.component_list)
                for i in b.params.component_list)
        blk._am_eq = Expression(blk.params.phase_list, rule=rule_am_eq)

        def rule_bm(b, p):
            return sum(b.mole_frac_phase_comp[p, i]*b.b[i]
                       for i in b.params.component_list)
        blk.bm = Expression(blk.params.phase_list, rule=rule_bm)

        def rule_A(b, p):
            return (b.am[p]*b.pressure /
                    (const.gas_constant*b.temperature)**2)
        blk.A = Expression(blk.params.phase_list, rule=rule_A)

        def rule_B(b, p):
            return (b.bm[p]*b.pressure /
                    (const.gas_constant*b.temperature))
        blk.B = Expression(blk.params.phase_list, rule=rule_B)

        def rule_A_eq(b, p):
            return (b._am_eq[p]*b.pressure /
                    (const.gas_constant*b._teq)**2)
        blk._A_eq = Expression(blk.params.phase_list, rule=rule_A_eq)

def supcost_rule(m, k):
        return sum(m.cs * m.s[k, j, t] * TimeStep for j in m.SUP for t in m.TIME.get_finite_elements())

    model.supcost = aml.Expression(model.SCEN, rule=supcost_rule)

    def boostcost_rule(m, k):
        return sum(m.ce * m.pow[k, j, t] * TimeStep for j in m.LINK_A for t in m.TIME.get_finite_elements())

    model.boostcost = aml.Expression(model.SCEN, rule=boostcost_rule)

    def trackcost_rule(m, k):
        return sum(
            m.cd * (m.dem[k, j, t] - m.stochd[k, j, t]) ** 2.0 for j in m.DEM for t in m.TIME.get_finite_elements())

    model.trackcost = aml.Expression(model.SCEN, rule=trackcost_rule)

    def sspcost_rule(m, k):
        return sum(
            m.cT * (m.px[k, i, m.TIME.last(), j] - m.px[k, i, m.TIME.first(), j]) ** 2.0 for i in m.LINK for j in m.DIS)

    model.sspcost = aml.Expression(model.SCEN, rule=sspcost_rule)

    def ssfcost_rule(m, k):
        return sum(
            m.cT * (m.fx[k, i, m.TIME.last(), j] - m.fx[k, i, m.TIME.first(), j]) ** 2.0 for i in m.LINK for j in m.DIS)

    model.ssfcost = aml.Expression(model.SCEN, rule=ssfcost_rule)

    def cost_rule(m, k):
        return 1e-6 * (m.supcost[k] + m.boostcost[k] + m.trackcost[k] + m.sspcost[k] + m.ssfcost[k])

blk.bm = Expression(blk.params.phase_list, rule=rule_bm)

        def rule_A(b, p):
            return (b.am[p]*b.pressure /
                    (const.gas_constant*b.temperature)**2)
        blk.A = Expression(blk.params.phase_list, rule=rule_A)

        def rule_B(b, p):
            return (b.bm[p]*b.pressure /
                    (const.gas_constant*b.temperature))
        blk.B = Expression(blk.params.phase_list, rule=rule_B)

        def rule_A_eq(b, p):
            return (b._am_eq[p]*b.pressure /
                    (const.gas_constant*b._teq)**2)
        blk._A_eq = Expression(blk.params.phase_list, rule=rule_A_eq)

        def rule_B_eq(b, p):
            return (b.bm[p]*b.pressure /
                    (const.gas_constant*b._teq))
        blk._B_eq = Expression(blk.params.phase_list, rule=rule_B_eq)

        blk.proc_Z_liq = ExternalFunction(library=_so,
                                          function="ceos_z_liq")
        blk.proc_Z_vap = ExternalFunction(library=_so,
                                          function="ceos_z_vap")
        blk.proc_Z_liq_x = ExternalFunction(library=_so,
                                            function="ceos_z_liq_extend")
        blk.proc_Z_vap_x = ExternalFunction(library=_so,
                                            function="ceos_z_vap_extend")

        def rule_delta(b, p, i):

def _enth_mol_phase(self):
        try:
            def rule_enth_mol_phase(b, p):
                p_config = b.params.get_phase(p).config
                return p_config.equation_of_state.enth_mol_phase(b, p)
            self.enth_mol_phase = Expression(self.params.phase_list,
                                             rule=rule_enth_mol_phase)
        except AttributeError:
            self.del_component(self.enth_mol_phase)
            raise

doc='EoS S factor')

        def func_b(b, j):
            return b.EoS_Bc*const.gas_constant *\
                   b.params.temperature_crit[j]/b.params.pressure_crit[j]
        blk.b = Param(blk.params.component_list,
                      initialize=func_b,
                      doc='Component b coefficient')

        def func_a(b, j):
            return (b.omegaA*((const.gas_constant *
                               b.params.temperature_crit[j])**2 /
                              b.params.pressure_crit[j]) *
                    ((1+b.fw[j]*(1-sqrt(b.temperature /
                                        b.params.temperature_crit[j])))**2))
        blk.a = Expression(blk.params.component_list,
                           rule=func_a,
                           doc='Component a coefficient')

        def func_a_eq(b, j):
            return (b.omegaA*((const.gas_constant *
                               b.params.temperature_crit[j])**2 /
                              b.params.pressure_crit[j]) *
                    ((1+b.fw[j]*(1-sqrt(b._teq /
                                        b.params.temperature_crit[j])))**2))
        blk._a_eq = Expression(blk.params.component_list,
                               rule=func_a_eq,
                               doc='Component a coefficient at Teq')

        def rule_am(b, p):
            return sum(sum(
                b.mole_frac_phase_comp[p, i]*b.mole_frac_phase_comp[p, j] *

self.scaling_factor[self.phase_frac[p]] = 10
            if p == "Vap":
                return vf
            elif p == "Liq":
                return 1.0 - vf

        self.phase_frac = Expression(
            phlist, rule=rule_phase_frac, doc="Phase fraction [unitless]"
        )

        # Component flow (for units that need it)
        def component_flow(b, i):
            self.scaling_factor[self.flow_mol_comp[i]] = 1e-3
            return self.flow_mol

        self.flow_mol_comp = Expression(
            component_list,
            rule=component_flow,
            doc="Total flow (both phases) of component [mol/s]",
        )

        # Total (mixed phase) properties

        # Enthalpy
        if self.state_vars == StateVars.TPX:
            self.enth_mol = Expression(
                expr=sum(self.phase_frac[p] * self.enth_mol_phase[p] for p in phlist)
            )
            self.scaling_factor[self.enth_mol] = 1e-3
        # Internal Energy
        self.energy_internal_mol = Expression(
            expr=sum(

def _entr_mol_phase_comp(self):
        try:
            def rule_entr_mol_phase_comp(b, p, j):
                p_config = b.params.get_phase(p).config
                return p_config.equation_of_state.entr_mol_phase_comp(b, p, j)
            self.entr_mol_phase_comp = Expression(
                self.params._phase_component_set,
                rule=rule_entr_mol_phase_comp)
        except AttributeError:
            self.del_component(self.entr_mol_phase_comp)
            raise

def _mw_phase(self):
        try:
            def rule_mw_phase(b, p):
                return sum(b.mole_frac_phase_comp[p, j] *
                           b.params.get_component(j).mw_comp
                           for j in b.params.component_list)
            self.mw_phase = Expression(
                    self.params.phase_list,
                    doc="Average molecular weight of each phase",
                    rule=rule_mw_phase)
        except AttributeError:
            self.del_component(self.mw_phase)
            raise