<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 def add_thermodynamic(self, em=1000):
 """Apply thermodynamic constraints to the model. Adding these constraints restricts the solution space to only contain solutions that have no internal loops [Schilling00]_. This is solved as a MILP problem as described in [Muller13]_. The time to solve a problem with thermodynamic constraints is usually much longer than a normal FBA problem. The ``em`` parameter is the upper bound on the delta mu reaction variables. This parameter has to be balanced based on the model size since setting the value too low can result in the correct solutions being infeasible and setting the value too high can result in numerical instability which again makes the correct solutions infeasible. The default value should work in all cases as long as the model is not unusually large. """


        internal = set(r for r in self._model.reactions
                       if not self._model.is_exchange(r))

        # Reaction fluxes
        v = self._v

        # Indicator variable
        alpha = self._prob.namespace(internal, types=lp.VariableType.Binary)

        # Delta mu is the stoichiometrically weighted sum of the compound mus.
        dmu = self._prob.namespace(internal)

        for reaction_id in self._model.reactions:
            if not self._model.is_exchange(reaction_id):
                flux = v(reaction_id)
                alpha_r = alpha(reaction_id)
                dmu_r = dmu(reaction_id)

                lower, upper = self._model.limits[reaction_id]

                # Constrain the reaction to a direction determined by alpha
                # and contrain the delta mu to a value in [-em; -1] if
                # alpha is one, otherwise in [1; em].
                self._prob.add_linear_constraints(
                    flux >= lower * (1 - alpha_r),
                    flux <= upper * alpha_r,
                    dmu_r >= -em * alpha_r + (1 - alpha_r),
                    dmu_r <= em * (1 - alpha_r) - alpha_r)

        # Define mu variables
        mu = self._prob.namespace(self._model.compounds)

        tdbalance_lhs = {reaction_id: 0
                         for reaction_id in self._model.reactions}
        for spec, value in iteritems(self._model.matrix):
            compound, reaction_id = spec
            if not self._model.is_exchange(reaction_id):
                tdbalance_lhs[reaction_id] += mu(compound) * value
        for reaction_id, lhs in iteritems(tdbalance_lhs):
            if not self._model.is_exchange(reaction_id):
                self._prob.add_linear_constraints(lhs == dmu(reaction_id))