docs: 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))