diff --git a/Bioindustrial-Park b/Bioindustrial-Park
index ab0aa3ac..978706e6 160000
--- a/Bioindustrial-Park
+++ b/Bioindustrial-Park
@@ -1 +1 @@
-Subproject commit ab0aa3ac6d3b8510b230d7b9f0fca6479316843e
+Subproject commit 978706e698cdd88626d916a120c910066bc307e7
diff --git a/benchmark/acetic_acid_reactive_purification_system.py b/benchmark/acetic_acid_reactive_purification_system.py
new file mode 100644
index 00000000..59c0a81c
--- /dev/null
+++ b/benchmark/acetic_acid_reactive_purification_system.py
@@ -0,0 +1,178 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Sat Mar 16 13:38:11 2024
+
+@author: cortespea
+"""
+import thermosteam as tmo
+import biosteam as bst
+import numpy as np
+from thermosteam.constants import R
+from math import exp
+try:
+    from .profile import register
+except:
+    def register(*args, **kwargs):
+        return lambda f: f
+
+__all__ = (
+    'create_system_acetic_acid_reactive_purification',
+)
+
+@register(
+    'acetic_acid_reactive_purification', 'Acetic acid reactive purification',
+    10, [2, 4, 6, 8, 10], 'AA\nr. sep.'
+)
+def create_system_acetic_acid_reactive_purification(alg='sequential modular'):
+    bst.settings.set_thermo(['Water', 'AceticAcid', 'MethylAcetate', 'Methanol'], cache=True)
+    Gamma = tmo.equilibrium.DortmundActivityCoefficients(bst.settings.chemicals)
+    class Esterification(bst.KineticReaction):
+        def volume(self, stream): # kg of catalyst
+            rho_cat = 770 # kg / m3
+            liquid_volume = self.liquid_volume
+            catalyst_volume = 0.5 * liquid_volume
+            catalyst_mass = catalyst_volume * rho_cat
+            return catalyst_mass
+        
+        def rate(self, stream): # kmol/kg-catalyst/hr
+            K_AcOH = 3.15
+            K_MeOH = 5.64
+            K_MeOAc = 4.15
+            K_H2O = 5.24
+            K_adsorption = np.array([K_H2O, K_AcOH, K_MeOAc, K_MeOH])
+            MWs = stream.chemicals.MW
+            k0f = 8.497e6
+            EAf = 60.47
+            k0r = 6.127e5
+            EAr = 63.73
+            T = stream.T
+            mol = stream.mol.to_array()
+            F_mol = stream.F_mol
+            if not F_mol: return 0
+            x = mol / F_mol
+            gamma = Gamma(x, T)
+            activity = gamma * x
+            a = K_adsorption * activity / MWs
+            H2O, LA, BuLA, BuOH = activity
+            kf = k0f * exp(-EAf / (R * T))
+            kr = k0r * exp(-EAr / (R * T))
+            a_sum = sum(a)
+            return 3600 * (kf * LA * BuOH - kr * BuLA * H2O) / (a_sum * a_sum) # kmol / kg-catalyst / hr
+    
+    with bst.System(algorithm=alg) as sys:
+        feed = bst.Stream(
+            'feed',
+            LacticAcid=4.174,
+            Water=5.470,
+        )
+        feed.vle(V=0.5, P=101325)
+        makeup_butanol = bst.Stream('makeup_butanol', Butanol=1)
+        makeup_butanol.T = makeup_butanol.bubble_point_at_P().T
+        recycle_butanol = bst.Stream('recycle_butanol')
+        esterification_reflux = bst.Stream('esterification_reflux')
+        esterification = bst.MESHDistillation(
+            'esterification',
+            ins=(feed, makeup_butanol, recycle_butanol, esterification_reflux), 
+            outs=('empty', 'bottoms', 'esterification_distillate'),
+            N_stages=6,
+            feed_stages=(3, 3, 3, 0),
+            stage_specifications={
+                5: ('Boilup', 1),
+                0: ('Boilup', 0),
+            },
+            liquid_side_draws={
+                0: 1,
+            },
+            stage_reactions={
+                i: Esterification('LacticAcid + Butanol -> Water + ButylLactate', reactant='LacticAcid')
+                for i in range(1, 5)
+            },
+            maxiter=10,
+            # LHK=('LacticAcid', 'Butanol'),
+            LHK=('Butanol', 'ButylLactate'),
+            P=0.3 * 101325,
+        )
+        @esterification.add_specification(run=True)
+        def adjust_flow():
+            target = 5.85
+            makeup_butanol.imol['Butanol'] = max(target - recycle_butanol.imol['Butanol'], 0)
+        
+        esterification.simulate()
+        for i in esterification.stages: print(i.Hnet) 
+        esterification.show()
+        breakpoint()
+        # esterification.simulate()
+        # for i in esterification.stages: print(i.Hnet) 
+        # esterification.stage_reactions={
+        #         i: Esterification('LacticAcid + Butanol -> Water + ButylLactate', reactant='LacticAcid')
+        #         for i in range(1, 17)
+        #     }
+        # esterification.LHK=('Butanol', 'ButylLactate')
+        # breakpoint()
+        # esterification.simulate()
+        esterification_settler = bst.StageEquilibrium(
+            'esterification_settler',
+            ins=(esterification-2), 
+            outs=(esterification_reflux, 'water_rich'),
+            phases=('L', 'l'),
+            top_chemical='Butanol',
+        )
+        water_distiller = bst.BinaryDistillation(
+            ins=esterification_settler-1, outs=('water_rich_azeotrope', 'water'),
+            x_bot=0.0001, y_top=0.2, k=1.2, Rmin=0.01,
+            LHK=('Butanol', 'Water'),
+        )
+        splitter = bst.Splitter(ins=water_distiller-1, split=0.5) # TODO: optimize split
+        hydrolysis_reflux = bst.Stream('hydrolysis_reflux')
+        hydrolysis = bst.MESHDistillation(
+            'hydrolysis',
+            ins=(esterification-1, splitter-0, hydrolysis_reflux),
+            outs=('empty', 'lactic_acid', 'hydrolysis_distillate'),
+            N_stages=53,
+            feed_stages=(27, 50, 0),
+            stage_specifications={
+                0: ('Boilup', 0),
+                52: ('Boilup', 1),
+            },
+            liquid_side_draws={
+                0: 1.0,
+            },
+            stage_reactions={
+                i: Esterification('LacticAcid + Butanol -> Water + ButylLactate', reactant='LacticAcid')
+                for i in range(1, 52) # It will run in reverse
+            },
+            P=101325,
+            LHK=('Butanol', 'LacticAcid'),
+        )
+        
+        # @esterification.add_specification(run=True)
+        # def adjust_flow():
+        #     target = 5.85
+        #     makeup_butanol.imol['Butanol'] = max(target - recycle_butanol.imol['Butanol'], 0)
+        
+        # Decanter
+        butanol_rich_azeotrope = bst.Stream('butanol_rich_azeotrope')
+        hydrolysis_settler = bst.StageEquilibrium(
+            'settler',
+            ins=(hydrolysis-2, water_distiller-0, butanol_rich_azeotrope), 
+            outs=('butanol_rich_extract', hydrolysis_reflux),
+            phases=('L', 'l'),
+            top_chemical='Butanol',
+            T=310,
+        )
+        
+        # Butanol purification
+        butanol_distiller = bst.BinaryDistillation(
+            ins=(hydrolysis_settler-0),
+            outs=(butanol_rich_azeotrope, recycle_butanol),
+            x_bot=0.0001, y_top=0.6, k=1.2, Rmin=0.01,
+            LHK=('Water', 'Butanol'),
+        )
+        
+    return sys
+
+if __name__ == '__main__':
+    sys = create_system_acetic_acid_reactive_purification()
+    sys.flatten()
+    sys.diagram()
+    sys.simulate()
\ No newline at end of file
diff --git a/benchmark/lactic_acid_purification_system.py b/benchmark/lactic_acid_purification_system.py
index 5c2f9837..09e76edf 100644
--- a/benchmark/lactic_acid_purification_system.py
+++ b/benchmark/lactic_acid_purification_system.py
@@ -39,7 +39,7 @@ def rate(self, stream): # kmol/kg-catalyst/hr
             kf = 2.59e4 * exp(-5.340e4 / (R * T))
             kr = 3.80e3 * exp(-5.224e4 / (R * T))
             H2O, LA, BuLA, BuOH = stream.mol / stream.F_mol
-            return self.stoichiometry * 3600 * (kf * LA * BuOH - kr * BuLA * H2O) # kmol / kg-catalyst / hr
+            return 3600 * (kf * LA * BuOH - kr * BuLA * H2O) # kmol / kg-catalyst / hr
     
     with bst.System(algorithm=alg) as sys:
         feed = bst.Stream(
@@ -47,27 +47,30 @@ def rate(self, stream): # kmol/kg-catalyst/hr
             LacticAcid=4.174,
             Water=5.470,
         )
-        makeup_butanol = bst.Stream('makeup_butanol')
+        feed.vle(V=0.5, P=101325)
+        makeup_butanol = bst.Stream('makeup_butanol', Butanol=1)
+        makeup_butanol.T = makeup_butanol.bubble_point_at_P().T
         recycle_butanol = bst.Stream('recycle_butanol')
         esterification_reflux = bst.Stream('esterification_reflux')
         esterification = bst.MESHDistillation(
             'esterification',
             ins=(feed, makeup_butanol, recycle_butanol, esterification_reflux), 
             outs=('empty', 'bottoms', 'esterification_distillate'),
-            N_stages=17,
-            feed_stages=(1, 16, 16, 0),
+            N_stages=6,
+            feed_stages=(3, 3, 3, 0),
             stage_specifications={
-                16: ('Boilup', 1),
+                5: ('Boilup', 1),
                 0: ('Boilup', 0),
             },
             liquid_side_draws={
-                0: 1.0,
+                0: 1,
             },
             stage_reactions={
                 i: Esterification('LacticAcid + Butanol -> Water + ButylLactate', reactant='LacticAcid')
-                for i in range(1, 17)
+                for i in range(1, 5)
             },
-            maxiter=200,
+            maxiter=10,
+            # LHK=('LacticAcid', 'Butanol'),
             LHK=('Butanol', 'ButylLactate'),
             P=0.3 * 101325,
         )
@@ -75,7 +78,20 @@ def rate(self, stream): # kmol/kg-catalyst/hr
         def adjust_flow():
             target = 5.85
             makeup_butanol.imol['Butanol'] = max(target - recycle_butanol.imol['Butanol'], 0)
-            
+        
+        esterification.simulate()
+        for i in esterification.stages: print(i.Hnet) 
+        esterification.show()
+        breakpoint()
+        # esterification.simulate()
+        # for i in esterification.stages: print(i.Hnet) 
+        # esterification.stage_reactions={
+        #         i: Esterification('LacticAcid + Butanol -> Water + ButylLactate', reactant='LacticAcid')
+        #         for i in range(1, 17)
+        #     }
+        # esterification.LHK=('Butanol', 'ButylLactate')
+        # breakpoint()
+        # esterification.simulate()
         esterification_settler = bst.StageEquilibrium(
             'esterification_settler',
             ins=(esterification-2), 
diff --git a/biosteam/_settings.py b/biosteam/_settings.py
index bcffdba7..a4a03d42 100644
--- a/biosteam/_settings.py
+++ b/biosteam/_settings.py
@@ -32,7 +32,7 @@ def {flowmethname}(self):
 get_unit_utility_cost_executable = '''
 def {costmethname}(self):
     """Return the {docname} cost [USD/hr]."""
-    return bst.stream_prices[name] * self.{flowmethname}()
+    return bst.stream_prices[{name}] * self.{flowmethname}()
 
 bst.Unit.{costmethname} = {costmethname}
 '''
@@ -199,7 +199,8 @@ def define_allocation_property(
 
 # %% Register stream utilities
 
-settings.register_utility('Natural gas', 0.218)
+settings.register_utility('Fuel', 0.218) # Natural gas
+settings.register_utility('Natural gas', 0.218) # Natural gas
 settings.register_utility('Ash disposal', -0.0318)
 settings.register_utility('Reverse osmosis water', 5.6e-4)
 settings.register_utility('Process water', 2.7e-4)
diff --git a/biosteam/facilities/_boiler_turbogenerator.py b/biosteam/facilities/_boiler_turbogenerator.py
index ae8af6d4..9b2e59da 100644
--- a/biosteam/facilities/_boiler_turbogenerator.py
+++ b/biosteam/facilities/_boiler_turbogenerator.py
@@ -40,7 +40,7 @@ class BoilerTurbogenerator(bst.Facility):
         * [0] Liquid/solid feed that will be burned.
         * [1] Gas feed that will be burned.
         * [2] Make-up water. 
-        * [3] Natural gas to satisfy steam and power requirement.
+        * [3] Natural gas/fuel to satisfy steam and electricity demand.
         * [4] Lime for flue gas desulfurization.
         * [5] Boiler chemicals.
     outs : 
@@ -55,7 +55,9 @@ class BoilerTurbogenerator(bst.Facility):
         Steam produced. Defaults to low pressure steam.
     other_agents = () : Iterable[UtilityAgent], optional
         Other steams produced. Defaults to all other heating agents.
-    natural_gas_price : float, optional
+    fuel_source : str, optional
+        Name fuel used to satisfy steam and electricity demand. Defaults to 'CH4'.
+    fuel_price : float, optional
         Price of natural gas [USD/kg]. Same as `bst.stream_utility_prices['Natural gas']`,
         which defaults to 0.218.
     ash_disposal_price : float, optional
@@ -159,21 +161,25 @@ def __init__(self, ID='', ins=None,
                  side_steam=None,
                  agent=None,
                  other_agents=None,
+                 fuel_price=None,
                  natural_gas_price=None,
                  ash_disposal_price=None,
                  T_emissions=None,
                  satisfy_system_electricity_demand=None,
                  boiler_efficiency_basis=None,
+                 fuel_source=None,
         ):
         if boiler_efficiency_basis is None: boiler_efficiency_basis = 'LHV'
         if boiler_efficiency is None: boiler_efficiency = 0.80
         if turbogenerator_efficiency is None: turbogenerator_efficiency = 0.85
         if satisfy_system_electricity_demand is None: satisfy_system_electricity_demand = True
+        if fuel_source is None: fuel_source = 'CH4'
         bst.Facility.__init__(self, ID, ins, outs, thermo)
         settings = bst.settings
         self.boiler_efficiency_basis = boiler_efficiency_basis
         self.agent = agent = agent or settings.get_heating_agent('low_pressure_steam')
-        self.define_utility('Natural gas', self.natural_gas)
+        self.fuel_source = fuel_source
+        self.define_utility('Fuel', self.fuel)
         self.define_utility('Ash disposal', self.ash_disposal)
         self.boiler_efficiency = boiler_efficiency
         self.turbogenerator_efficiency = turbogenerator_efficiency
@@ -183,7 +189,10 @@ def __init__(self, ID='', ins=None,
         self.side_steam = side_steam
         self.other_agents = [i for i in settings.heating_agents if i is not agent] if other_agents is None else other_agents
         self.T_emissions = self.agent.T if T_emissions is None else T_emissions # Assume no heat integration
-        if natural_gas_price is not None: self.natural_gas_price = natural_gas_price
+        if natural_gas_price is not None: 
+            self.fuel_price = natural_gas_price
+        elif fuel_price is not None: 
+            self.fuel_price = fuel_price
         if ash_disposal_price is not None: self.ash_disposal_price = ash_disposal_price
         self.satisfy_system_electricity_demand = satisfy_system_electricity_demand
       
@@ -226,9 +235,10 @@ def makeup_water(self):
         return self.ins[2]
     
     @property
-    def natural_gas(self):
-        """[Stream] Natural gas to satisfy steam and electricity requirements."""
+    def fuel(self):
+        """[Stream] Fuel used to satisfy steam and electricity requirements."""
         return self.ins[3]
+    natural_gas = fuel
     
     @property
     def ash_disposal(self):
@@ -236,13 +246,13 @@ def ash_disposal(self):
         return self.outs[2]
     
     @property
-    def natural_gas_price(self):
+    def fuel_price(self):
         """[Float] Price of natural gas, same as `bst.stream_utility_prices['Natural gas']`."""
-        return bst.stream_utility_prices['Natural gas']
-    
-    @natural_gas_price.setter
-    def natural_gas_price(self, new_price):
-        bst.stream_utility_prices['Natural gas'] = new_price
+        return bst.stream_utility_prices['Fuel']
+    @fuel_price.setter
+    def fuel_price(self, new_price):
+        bst.stream_utility_prices['Fuel'] = new_price
+    natural_gas_price = fuel_price
     
     @property
     def ash_disposal_price(self):
@@ -277,10 +287,11 @@ def _design(self):
         chemicals = self.chemicals
         self._load_utility_agents()
         mol_steam = sum([i.flow for i in self.steam_utilities])
-        feed_solids, feed_gas, makeup_water, feed_CH4, lime, chems = self.ins
-        feed_CH4.phase = 'g'
-        feed_CH4.set_property('T', 60, 'degF')
-        feed_CH4.set_property('P', 14.73, 'psi')
+        feed_solids, feed_gas, makeup_water, fuel, lime, chems = self.ins
+        if self.fuel_source == 'CH4':
+            fuel.phase = 'g'
+            fuel.set_property('T', 60, 'degF')
+            fuel.set_property('P', 14.73, 'psi')
         emissions, blowdown_water, ash_disposal = self.outs
         if not lime.price:
             lime.price = 0.19937504680689402
@@ -300,21 +311,22 @@ def _design(self):
         self.combustion_reactions = combustion_rxns = chemicals.get_combustion_reactions()
         non_empty_feeds = [i for i in (feed_solids, feed_gas) if not i.isempty()]
         boiler_efficiency_basis = self.boiler_efficiency_basis
-        def calculate_excess_electricity_at_natual_gas_flow(natural_gas_flow):
-            if natural_gas_flow:
-                natural_gas_flow = abs(natural_gas_flow)
-                feed_CH4.imol['CH4'] = natural_gas_flow
+        fuel_source = self.fuel_source
+        def calculate_excess_electricity_at_natual_gas_flow(fuel_flow):
+            if fuel_flow:
+                fuel_flow = abs(fuel_flow)
+                fuel.imol[fuel_source] = fuel_flow
             else:
-                feed_CH4.empty()
-            emissions_mol[:] = feed_CH4.mol
+                fuel.empty()
+            emissions_mol[:] = fuel.mol
             for feed in non_empty_feeds: emissions_mol[:] += feed.mol
             combustion_rxns.force_reaction(emissions_mol)
             emissions.imol['O2'] = 0
             if boiler_efficiency_basis == 'LHV':
-                H_combustion = feed_CH4.LHV
+                H_combustion = fuel.LHV
                 for feed in non_empty_feeds: H_combustion += feed.LHV
             elif boiler_efficiency_basis == 'HHV':
-                H_combustion = feed_CH4.HHV
+                H_combustion = fuel.HHV
                 for feed in non_empty_feeds: H_combustion += feed.HHV
             else:
                 raise ValueError(
@@ -339,11 +351,12 @@ def calculate_excess_electricity_at_natual_gas_flow(natural_gas_flow):
             else:
                 return work
         
-        self._excess_electricity_without_natural_gas = excess_electricity = calculate_excess_electricity_at_natual_gas_flow(0)
+        self._excess_electricity_without_fuel = excess_electricity = calculate_excess_electricity_at_natual_gas_flow(0)
         if excess_electricity < 0:
             f = calculate_excess_electricity_at_natual_gas_flow
             lb = 0.
-            ub = - excess_electricity * 3600 / feed_CH4.chemicals.CH4.LHV
+            fuel.imol[fuel_source] = 1
+            ub = - excess_electricity * 3600 / fuel.LHV
             while f(ub) < 0.: 
                 lb = ub
                 ub *= 2
diff --git a/biosteam/plots/plots.py b/biosteam/plots/plots.py
index abb42d21..b18fd99e 100644
--- a/biosteam/plots/plots.py
+++ b/biosteam/plots/plots.py
@@ -699,15 +699,17 @@ def plot_spearman_2d(rhos, top=None, name=None, color_wheel=None, index=None,
 class Box:
     __slots__ = ('axis', 'light', 'dark', '_position', '_baseline_position')
     
-    def __init__(self, axis, position=0, light=None, dark=None):
+    def __init__(self, axis, position=None, light=None, dark=None):
         self.axis = axis
         self.light = light
         self.dark = dark
-        self._baseline_position = self._position = position
+        if position is None: position = [0]
+        self._baseline_position = position[0]
+        self._position = position
         
     def get_position(self, shift=1):
-        self._position += shift
-        return self._position
+        self._position[0] += shift
+        return self._position[0]
     
     def reset(self):
         self._position = self._baseline_position
@@ -849,6 +851,10 @@ def plot_kde(x, y, nbins=100, ax=None, fig=None,
              aspect_ratio=1.25, cmaps=None, xbox_width=None,
              ybox_width=None, **kwargs):
     axis_not_given = ax is None
+    xs = x if isinstance(x, (tuple, list)) else (x,)
+    ys = y if isinstance(y, (tuple, list)) else (y,)
+    N_xs = len(xs)
+    N_ys = len(ys)
     if axis_not_given:
         grid_kw = dict(height_ratios=[1, 8], width_ratios=[8, aspect_ratio])
         fig, all_axes = plt.subplots(
@@ -859,12 +865,30 @@ def plot_kde(x, y, nbins=100, ax=None, fig=None,
         ax = all_axes[1, 0]
         xbox_ax = all_axes[0, 0]
         ybox_ax = all_axes[1, 1]
-        xbox = Box(xbox_ax, **(xbox_kwargs or {}))
-        ybox = Box(ybox_ax, **(ybox_kwargs or {}))
-    xs = x if isinstance(x, tuple) else (x,)
-    ys = y if isinstance(y, tuple) else (y,)
+    if xbox is None: 
+        if xbox_kwargs is None: xbox_kwargs = {}
+        position = [0]
+        if isinstance(xbox_kwargs, dict):
+            xboxes = [Box(xbox_ax, position, **xbox_kwargs)] * N_xs
+        else:
+            xboxes = [Box(xbox_ax, position, **i) for i in xbox_kwargs]
+    elif isinstance(xbox, Box):
+        xboxes = [xbox] * N_xs
+    else:
+        xboxes = xbox
+    if ybox is None:
+        if ybox_kwargs is None: ybox_kwargs = {}
+        position = [0]
+        if isinstance(ybox_kwargs, dict):
+            yboxes = [Box(ybox_ax, position, **ybox_kwargs)] * N_xs
+        else:
+            yboxes = [Box(ybox_ax, position, **i) for i in ybox_kwargs]
+    elif isinstance(ybox, Box):
+        yboxes = [ybox] * N_ys
+    else:
+        yboxes = ybox
     if cmaps is None: cmaps = len(xs) * [None]
-    for x, y, cmap in zip(xs, ys, cmaps):
+    for x, y, cmap, xbox, ybox in zip(xs, ys, cmaps, xboxes, yboxes):
         # Evaluate a gaussian kde on a regular grid of nbins x nbins over data extents
         k = kde.gaussian_kde([x, y])
         z = k(np.vstack([x, y]))
diff --git a/biosteam/plots/utils.py b/biosteam/plots/utils.py
index cb46cccb..5bdfa5d5 100644
--- a/biosteam/plots/utils.py
+++ b/biosteam/plots/utils.py
@@ -75,9 +75,9 @@ def colorbar(self, fig, ax, colorplot, label=True, **cbarkwargs):
         # cbar_ax.locator_params(nbins=self.N_ticks)
         if label:
             if self.title_position == 'top':
-                cbar_ax.set_title(self.title, self.ylabelkwargs)
+                cbar_ax.set_title(self.title, **self.ylabelkwargs)
             else:
-                cbar_ax.set_ylabel(self.title, self.ylabelkwargs)
+                cbar_ax.set_ylabel(self.title, **self.ylabelkwargs)
         try:
             ylabels = [y.get_text() for y in cbar_ax.get_yticklabels()]
             ylabels = [(i if i[0].isdigit() else '-'+i[1:]) for i in ylabels]
diff --git a/biosteam/process_tools/process_model.py b/biosteam/process_tools/process_model.py
index aa8a41f4..5d1cfaac 100644
--- a/biosteam/process_tools/process_model.py
+++ b/biosteam/process_tools/process_model.py
@@ -46,4 +46,5 @@ def _info(self, p, m):
     def show(self, p=None, m=None):
         """Return information on p-parameters and m-metrics."""
         print(self._info(p, m))
-    _ipython_display_ = show
\ No newline at end of file
+    _ipython_display_ = show
+    
\ No newline at end of file
diff --git a/biosteam/units/aerated_bioreactor.py b/biosteam/units/aerated_bioreactor.py
index fd0c4d7e..3139f501 100644
--- a/biosteam/units/aerated_bioreactor.py
+++ b/biosteam/units/aerated_bioreactor.py
@@ -21,7 +21,7 @@
 from scipy.constants import g
 import flexsolve as flx
 from warnings import filterwarnings, catch_warnings
-from scipy.optimize import minimize_scalar, minimize, least_squares
+from scipy.optimize import minimize_scalar, minimize, least_squares, differential_evolution
 from biosteam.units.design_tools import aeration
 
 __all__ = (
@@ -90,7 +90,7 @@ class AeratedBioreactor(StirredTankReactor):
     batch_default = True
     default_methods = {
         'Stirred tank': 'Riet',
-        'Bubble column': 'DeJesus',
+        'Bubble column': 'Dewes',
     }
     def _init(
             self, reactions, theta_O2=0.5, Q_O2_consumption=None,
@@ -129,19 +129,31 @@ def get_kLa(self):
             operating_time = self.tau / self.design_results.get('Batch time', 1.)
             N_reactors = self.parallel['self']
             P = 1000 * self.kW_per_m3 * V # W
-            air_in = self.cooled_compressed_air
+            air_in = self.sparged_gas
             D = self.get_design_result('Diameter', 'm')
             F = air_in.get_total_flow('m3/s') / N_reactors / operating_time
             R = 0.5 * D
             A = pi * R * R
             self.superficial_gas_flow = U = F / A # m / s 
             return aeration.kLa_stirred_Riet(P, V, U, **self.kLa_kwargs) # 1 / s 
+        elif self.kLa is aeration.kla_bubcol_Dewes:
+            V = self.get_design_result('Reactor volume', 'm3') * self.V_wf
+            operating_time = self.tau / self.design_results.get('Batch time', 1.)
+            N_reactors = self.parallel['self']
+            air_in = self.sparged_gas
+            D = self.get_design_result('Diameter', 'm')
+            F = air_in.get_total_flow('m3/s') / N_reactors / operating_time
+            R = 0.5 * D
+            A = pi * R * R
+            self.superficial_gas_flow = U = F / A # m / s 
+            feed = self.ins[0]
+            return aeration.kla_bubcol_Dewes(U, feed.get_property('mu', 'mPa*s'), air_in.get_property('rho', 'kg/m3'))
         else:
             raise NotImplementedError('kLa method has not been implemented in BioSTEAM yet')
     
     def get_agitation_power(self, kLa):
         if self.kLa is aeration.kLa_stirred_Riet:
-            air_in = self.cooled_compressed_air
+            air_in = self.sparged_gas
             N_reactors = self.parallel['self']
             operating_time = self.tau / self.design_results.get('Batch time', 1.)
             V = self.get_design_result('Reactor volume', 'm3') * self.V_wf
@@ -174,7 +186,7 @@ def air(self):
             if i.phase == 'g': return i
     
     @property
-    def cooled_compressed_air(self):
+    def sparged_gas(self):
         return self.air_cooler.outs[0]
     
     @property
@@ -216,7 +228,7 @@ def _run(self):
             if OUR > 0: effluent.imol['O2'] = 0.
             self._run_vent(vent, effluent)
             return
-        air_cc = self.cooled_compressed_air
+        air_cc = self.sparged_gas
         air_cc.copy_like(air)
         air_cc.P = compressor.P = self._inlet_air_pressure()
         air_cc.T = self.T
@@ -257,7 +269,7 @@ def _run_reactions(self, effluent):
         self.reactions.force_reaction(effluent)
     
     def _solve_total_power(self, OUR): # For OTR = OUR [mol / s]
-        air_in = self.cooled_compressed_air
+        air_in = self.sparged_gas
         N_reactors = self.parallel['self']
         operating_time = self.tau / self.design_results.get('Batch time', 1.)
         V = self.get_design_result('Reactor volume', 'm3') * self.V_wf
@@ -289,7 +301,7 @@ def get_OTR(self):
         operating_time = self.tau / self.design_results.get('Batch time', 1.)
         N_reactors = self.parallel['self']
         kLa = self.get_kLa()
-        air_in = self.cooled_compressed_air
+        air_in = self.sparged_gas
         vent = self.vent
         P_O2_air = air_in.get_property('P', 'bar') * air_in.imol['O2'] / air_in.F_mol
         P_O2_vent = vent.get_property('P', 'bar') * vent.imol['O2'] / vent.F_mol
@@ -435,7 +447,10 @@ def _init(self,
         self.mixins = {} if mixins is None else mixins # dict[int, tuple[int]] Pairs of variable feed gas index and inlets that will be mixed.
         StirredTankReactor._init(self, **kwargs)
         if design is None: 
-            design = 'Stirred tank'
+            if self.kW_per_m3 == 0:
+                design = 'Bubble column'
+            else:
+                design = 'Stirred tank'
         elif design not in aeration.kLa_method_names:
             raise ValueError(
                 f"{design!r} is not a valid design; only "
@@ -565,6 +580,7 @@ def load_flow_rates(F_feeds):
             run_auxiliaries()
             effluent.set_data(effluent_liquid_data)
             effluent.mix_from([self.sparged_gas, -s_consumed, s_produced, *liquid_feeds], energy_balance=False)
+            if (effluent.mol < 0).any(): breakpoint()
             vent.empty()
             self._run_vent(vent, effluent)
         
@@ -574,7 +590,7 @@ def load_flow_rates(F_feeds):
             breakpoint()
             bst.Stream.sum(self.normal_gas_feeds, energy_balance=False)
         baseline_flows = baseline_feed.get_flow('mol/s', self.gas_substrates)
-        bounds = np.array([[max(SURs[i] - baseline_flows[i], 0), 5 * SURs[i]] for i in index])
+        bounds = np.array([[max(1.01 * SURs[i] - baseline_flows[i], 0), 10 * SURs[i]] for i in index])
         if self.optimize_power:
             def total_power_at_substrate_flow(F_substrates):
                 load_flow_rates(F_substrates)
@@ -595,13 +611,20 @@ def gas_flow_rate_objective(F_substrates):
                 mask = STRs - F_ins > 0
                 STRs[mask] = F_ins[mask]
                 diff = SURs - STRs
+                diff[diff > 0] *= 1e6 # Force transfer rate to meet uptake rate
                 return (diff * diff).sum()
             
             f = gas_flow_rate_objective
-            bounds = bounds.T
             with catch_warnings():
                 filterwarnings('ignore')
+                bounds = bounds.T
                 results = least_squares(f, 1.2 * SURs, bounds=bounds, ftol=SURs.min() * 1e-6)
+                # results = differential_evolution(f, bounds=bounds, tol=SURs.min() * 1e-6)
+                # print('----')
+                # print(results)
+            if getattr(self, 'debug', None):
+                breakpoint()
+            self._results = results
             load_flow_rates(results.x / x_substrates)
             
         
@@ -641,14 +664,8 @@ def get_STRs(self):
         V = self.get_design_result('Reactor volume', 'm3') * self.V_wf
         operating_time = self.tau / self.design_results.get('Batch time', 1.)
         N_reactors = self.parallel['self']
-        P = 1000 * self.kW_per_m3 * V # W
         gas_in = self.sparged_gas
-        D = self.get_design_result('Diameter', 'm')
-        F = gas_in.get_total_flow('m3/s') / N_reactors / operating_time
-        R = 0.5 * D
-        A = pi * R * R
-        self.superficial_gas_flow = U = F / A # m / s 
-        kLa = aeration.kLa_stirred_Riet(P, V, U, **self.kLa_kwargs) # 1 / s 
+        kLa = self.get_kLa() # 1 / s 
         vent = self.vent
         P_gas = gas_in.get_property('P', 'bar')
         P_vent = vent.get_property('P', 'bar')
diff --git a/biosteam/units/design_tools/aeration.py b/biosteam/units/design_tools/aeration.py
index b7f7a410..41b6064a 100644
--- a/biosteam/units/design_tools/aeration.py
+++ b/biosteam/units/design_tools/aeration.py
@@ -10,6 +10,8 @@
 .. autofunction:: biosteam.units.design_tools.aeration.C_L
 .. autofunction:: biosteam.units.design_tools.aeration.C_O2_L
 .. autofunction:: biosteam.units.design_tools.aeration.kLa_stirred_Riet
+.. autofunction:: biosteam.units.design_tools.aeration.kla_stirred_Labik
+.. autofunction:: biosteam.units.design_tools.aeration.kLa_stirred_Galaction
 .. autofunction:: biosteam.units.design_tools.aeration.kla_bubcol_Deshpande
 .. autofunction:: biosteam.units.design_tools.aeration.kla_bubcol_DeJesus
 .. autofunction:: biosteam.units.design_tools.aeration.kla_bubcol_Akita_Yoshida
@@ -17,8 +19,7 @@
 .. autofunction:: biosteam.units.design_tools.aeration.kla_bubcol_Seno
 .. autofunction:: biosteam.units.design_tools.aeration.kla_bubcol_Suh
 .. autofunction:: biosteam.units.design_tools.aeration.kla_bubcol_Shah
-.. autofunction:: biosteam.units.design_tools.aeration.kla_stirred_Labik
-.. autofunction:: biosteam.units.design_tools.aeration.kLa_stirred_Galaction
+.. autofunction:: biosteam.units.design_tools.aeration.kla_bubcol_Dewes
 .. autofunction:: biosteam.units.design_tools.aeration.P_at_kLa_Riet
 
 References
@@ -56,6 +57,8 @@
     Prediction of oxygen mass transfer coefficients in stirred bioreactors for bacteria, yeasts and fungus broths.
     Biochemical Engineering Journal, 20(1), 85–94. https://doi.org/10.1016/j.bej.2004.02.005
 .. [11] Benz, G. T. Bioreactor Design for Chemical Engineers. AICHE CEP 2011, 21–26.
+.. [12] Dewes, I., & Schumpe, A. (1997). Gas density effect on mass transfer in the slurry bubble column. 
+    Chemical Engineering Science, 52(21–22), 4105–4109. https://doi.org/10.1016/S0009-2509(97)00252-2
 
 """
 
@@ -395,31 +398,30 @@ def kla_bubcol_Suh(D, mu_l, rho_l, D_l, g, sigma_l, epsilon_g, V_g, V_l, coeffic
 
     kla = (D_l/(D**2)) * 0.018 * (mu_l/(rho_l * D_l))**(0.5) * ((g*rho_l*D**2)/sigma_l)**(0.75) * (g*D**3*rho_l**2/(mu_l**2))**(0.39)*(V_g/(sqrt(g*D)))**0.51
     return kla
+
 @register
 def kla_bubcol_Dewes(V_g, mu_l, rho_g):
     """
-    Returns the KLa coefficient for a bubble column reactor based on the Dewes & Schumpe (1997) correlation.
-    -------
-    Parameters:
+    Returns the KLa coefficient for a bubble column reactor based on the Dewes & Schumpe (1997) correlation [12]_.
+    
+    Parameters
+    ----------
     V_g: float
         Superficial gas velocity, [m/s]
     mu_l: float
         Viscosity of the liquid, [mPa.s]
     rho_g: float
         Density of the gas, [kg/m^3]
-    -------
-    References: 
-    Dewes, I., & Schumpe, A. (1997). 
-        Gas density effect on mass transfer in the slurry bubble column. 
-        Chemical Engineering Science, 52(21–22), 4105–4109. https://doi.org/10.1016/S0009-2509(97)00252-2
-    -------
-    Notes:
+    
+    Notes
+    -----
     This correlation is valid under the following range of values:
     D_c: 0.115 m -> Diameter of the column [m]
     H_l: 1.37 m -> Liquid height in the column [m]
     V_g: 0.01 - 0.08 m/s -> Superficial gas velocity [m/s]
     mu_l: 1.35 - 583 mPa.s -> Viscosity of the liquid [mPa.s]
     rho_g: 0.4 - 18.8 kg/m^3 -> Density of the gas [kg/m^3]
+    
     """
 
     kla = V_g**0.90 * mu_l**(-0.55) * rho_g**0.46
diff --git a/biosteam/units/distillation.py b/biosteam/units/distillation.py
index f5f5a7c8..917782d3 100644
--- a/biosteam/units/distillation.py
+++ b/biosteam/units/distillation.py
@@ -2536,14 +2536,14 @@ def _init(self,
         self._last_args = (
             self.N_stages, self.feed_stages, self.vapor_side_draws, 
             self.liquid_side_draws, self.use_cache, *self._ins, 
-            self.partition_data, self.P,
-            reflux, boilup,
+            self.partition_data, self.P, self.stage_specifications,
         )
         
     @property
     def reflux(self):
-        name, value = self.stage_specifications[0]
-        if name == 'Reflux': return value
+        if 0 in self.stage_specifications:
+            name, value = self.stage_specifications[0]
+            if name == 'Reflux': return value
     @reflux.setter
     def reflux(self, reflux):
         self.stage_specifications[0] = ('Reflux', reflux)
@@ -2565,8 +2565,7 @@ def _setup(self):
         super()._setup()
         args = (self.N_stages, self.feed_stages, self.vapor_side_draws, 
                 self.liquid_side_draws, self.use_cache, *self._ins, 
-                self.partition_data, self.P, 
-                self.reflux, self.boilup)
+                self.partition_data, self.P, self.stage_specifications)
         if args != self._last_args:
             MultiStageEquilibrium._init(
                 self, N_stages=self.N_stages,
@@ -2668,8 +2667,8 @@ def _actual_stages(self):
             return self.N_stages, np.ceil(self.N_stages / eff)
        
     def update_liquid_holdup(self):
-        diameter = self.estimate_diameter()
-        hw = weir_height = diameter * self.weir_height * 12 # inches
+        diameter = self.estimate_diameter() 
+        hw = weir_height = self._TS * self.weir_height * 0.0393701 # mm to inches
         area = diameter * diameter * pi / 4
         b = 2/3
         partitions = self.partitions
@@ -2684,7 +2683,7 @@ def update_liquid_holdup(self):
             else:
                 rho_V = vapor.rho
                 rho_L = liquid.rho
-                active_area = area * (1 - partition.downcomer_area_fraction)
+                active_area = area * (1 - 2 * partition.downcomer_area_fraction)
                 Ua = vapor.get_total_flow('ft3/s') / active_area
                 Ks = Ua * sqrt(rho_V / (rho_L - rho_V)) # Capacity parameter    
             Phi_e = exp(-4.257 * Ks**0.91) # Effective relative froth density 
diff --git a/biosteam/units/stage.py b/biosteam/units/stage.py
index 67bf029f..17edb195 100644
--- a/biosteam/units/stage.py
+++ b/biosteam/units/stage.py
@@ -560,7 +560,7 @@ def _create_material_balance_equations(self):
     def _update_energy_variable(self, departure):
         phases = self.phases
         if phases == ('g', 'l'):
-            self.B += self.partition.B_relaxation_factor * departure
+            self.B += (1 - self.partition.B_relaxation_factor) * departure
         elif phases == ('L', 'l'):
             self.T = T = self.T + departure
             for i in self.outs: i.T = T
@@ -591,10 +591,10 @@ class PhasePartition(Unit):
     _N_ins = 1
     _N_outs = 2
     strict_infeasibility_check = False
-    dmol_relaxation_factor = 0.5
-    B_relaxation_factor = 0.5
-    T_relaxation_factor = 0.5
-    K_relaxation_factor = 0.5
+    dmol_relaxation_factor = 0.9
+    B_relaxation_factor = 0.9
+    T_relaxation_factor = 0.9
+    K_relaxation_factor = 0.9
     
     def _init(self, phases, partition_data, top_chemical=None, reaction=None):
         self.partition_data = partition_data
@@ -785,6 +785,11 @@ def _run_decoupled_reaction(self, P=None, relaxation_factor=None):
         f = self.dmol_relaxation_factor if relaxation_factor is None else relaxation_factor
         old = self.dmol
         new = self.reaction.conversion(bottom)
+        # mol = bottom.mol
+        # mask = (mol + old) < 0
+        # if mask.any(): 
+        #     old *= (-mol[mask] / old[mask]).min() * (1 - 1e-6)
+        #     assert (old + mol >= 0).all()
         self.dmol = f * old + (1 - f) * new
     
     def _run_lle(self, P=None, update=True, top_chemical=None):
@@ -1065,16 +1070,16 @@ class MultiStageEquilibrium(Unit):
     """
     _N_ins = 2
     _N_outs = 2
-    default_maxiter = 20
-    default_max_attempts = 20
-    default_fallback_maxiter = 1
+    max_attempts = 5
+    default_maxiter = 50
+    default_fallback = None
     default_molar_tolerance = 1e-3
     default_relative_molar_tolerance = 1e-6
     default_algorithm = 'root'
     available_algorithms = {'root', 'optimize'}
     default_methods = {
         'root': 'fixed-point',
-        'optimize': 'SLSQP',
+        'optimize': 'CG',
         'SurPASS': 'differential evolution', 
     }
     
@@ -1084,9 +1089,9 @@ class MultiStageEquilibrium(Unit):
         'wegstein': (flx.conditional_wegstein, True, {}),
     }
     optimize_options: dict[str, tuple[Callable, dict]] = {
-        'SLSQP': (minimize, True, True, {'tol': 0.1, 'method': 'SLSQP'}),
-        'CG': (minimize, False, False, {'tol': 1e-3, 'method': 'CG', 'options':{'maxiter': 500}}),
-        'differential evolution': (differential_evolution, False, True, {'seed': 0, 'popsize': 12, 'tol': 1e-6}),
+        'SLSQP': (minimize, True, True, {'tol': 0.1, 'method': 'SLSQP'}), # This one does not work
+        'CG': (minimize, False, False, {'tol': 1e-3, 'method': 'CG', 'options':{'maxiter': 500}}), # This one is great
+        'differential evolution': (differential_evolution, False, True, {'seed': 0, 'popsize': 12, 'tol': 1e-6}), # This works, but its extremely slow
     }
     SurPASS_options: dict[str, tuple[Callable, dict]] = {
         'differential evolution': (
@@ -1313,8 +1318,8 @@ def _init(self,
         #: [int] Maximum number of iterations.
         self.maxiter = self.default_maxiter if maxiter is None else maxiter
         
-        #: [int] Maximum number of iterations for fallback algorithm.
-        self.fallback_maxiter = self.default_fallback_maxiter
+        #: tuple[str, str] Fallback algorithm and method.
+        self.fallback = self.default_fallback
 
         #: [float] Molar tolerance (kmol/hr)
         self.molar_tolerance = self.default_molar_tolerance
@@ -1331,8 +1336,6 @@ def _init(self,
         self.method = self.default_methods[self.algorithm] if method is None else method
     
         self.inside_out = inside_out
-        
-        self.max_attempts = self.default_max_attempts
     
     @property
     def aggregated_stages(self):
@@ -1553,7 +1556,6 @@ def outlet_stages(self):
             return outlet_stages
     
     def correct_overall_mass_balance(self):
-        return
         outmol = sum([i.mol for i in self.outs])
         inmol = sum([i.mol for i in self.ins])
         stage_reactions = self.stage_reactions
@@ -1587,7 +1589,7 @@ def _feed_flows_and_conversion(self):
             dmol = p.dmol
             feed_mol = p.feed.mol
             pseudo_feed = feed_mol + dmol
-            mask = (pseudo_feed < 0) & (dmol < 0)
+            mask = (pseudo_feed < 0)
             if pseudo_feed[mask].any():
                 x = (-feed_mol[mask] / dmol[mask]).min()
                 assert 0 <= x <= 1
@@ -1627,46 +1629,6 @@ def set_flow_rates(self, top_flows):
                     partition.B = tnet / bnet
             for i in stage.splitters: i._run()
         
-    def set_flow_rates_old(self, top_flows):
-        top, bottom = self.multi_stream.phases
-        flow_tol = -1e-6 * self.multi_stream.mol
-        stages = self.stages
-        N_stages = self.N_stages
-        range_stages = range(N_stages)
-        index = self._update_index
-        top_flows[top_flows < 0.] = 0.
-        has_infeasible_flow = True
-        infeasible_checks = set()
-        while has_infeasible_flow:
-            has_infeasible_flow = False
-            for i in range_stages:
-                stage = stages[i]
-                partition = stage.partition
-                s_top, _ = partition.outs
-                s_top.mol[index] = top_flows[i]
-                if stage.top_split: stage.splitters[0]._run()
-            for i in range_stages:
-                stage = stages[i]
-                partition = stage.partition
-                s_top, s_bottom = partition.outs
-                bottom_flow = sum([i.mol for i in stage.ins]) - s_top.mol
-                mask = bottom_flow < 0.
-                if mask.any():
-                    has_infeasible_flow = (bottom_flow[mask] < flow_tol[mask]).any()
-                    if i not in infeasible_checks and has_infeasible_flow:
-                        infeasible_checks.add(i)
-                        infeasible_index, = np.where(mask[index])
-                        # TODO: Find algebraic solution to keeping top flow rates within feasible region.
-                        # This is only a temporary solution.
-                        infeasible_flow = bottom_flow[mask]
-                        top_flows[i, infeasible_index] += infeasible_flow
-                        break
-                    else:
-                        has_infeasible_flow = False
-                        bottom_flow[mask] = 0.
-                s_bottom.mol[:] = bottom_flow
-                if stage.bottom_split: stage.splitters[-1]._run()
-            
     def _run(self):
         if all([i.isempty() for i in self.ins]): 
             for i in self.outs: i.empty()
@@ -1674,63 +1636,41 @@ def _run(self):
         try:
             top_flow_rates = self.hot_start()
             algorithm = self.algorithm
+            if self._has_lle:
+                self._psuedo_equilibrium_options = tmo.equilibrium.LLE.pseudo_equilibrium_inner_loop_options.copy()
+                self._psuedo_equilibrium_options['xtol'] *= self.feed_flows.max()
             if algorithm == 'root':
-                self.converged = True
+                method = self.method
+                solver, conditional, options = self.root_options[method]
+                if not conditional: raise NotImplementedError(f'method {self.method!r} not implemented in BioSTEAM (yet)')
                 for i in range(self.max_attempts):
-                    self.attempt = i
                     self.iter = 0
-                    self.fallback_iter = 0
-                    method = self.method
-                    solver, conditional, options = self.root_options[method]
-                    if conditional:
-                        top_flow_rates = solver(self._conditional_iter, top_flow_rates)
-                    else:
-                        raise NotImplementedError(f'method {self.method!r} not implemented in BioSTEAM (yet)')
-                    if method == 'fixed-point' and self.iter == self.maxiter:
-                        top_flow_rates = flx.conditional_fixed_point(
-                            self._sequential_iter, 
-                            top_flow_rates,
-                        )
-                        if self.fallback_iter < self.fallback_maxiter:
-                            break
+                    top_flow_rates = solver(self._conditional_iter, top_flow_rates)
+                    if self.iter == self.maxiter:
+                        for i in self.stages: i._run()
+                        for i in reversed(self.stages): i._run()
                     else:
                         break
-                else:
-                    self.converged = False
+                if self.iter == self.maxiter and self.fallback and self.fallback != (self.algorithm, self.method):
+                    algorithm, method = self.algorithm, self.method
+                    self.algorithm, self.method = self.fallback
+                    try:
+                        self._run()
+                    finally:
+                        self.algorithm, self.method = algorithm, method
+                    self.iter = 0
+                    top_flow_rates = solver(self._conditional_iter, top_flow_rates)
+            elif algorithm == 'sequential':
+                self.iter = 0
+                top_flow_rates = flx.conditional_fixed_point(
+                    self._sequential_iter, 
+                    top_flow_rates,
+                )
             elif algorithm == 'optimize':
                 solver, constraints, bounded, options = self.optimize_options[self.method]
                 if constraints and bounded:
                     raise NotImplementedError(f'optimize method {self.method!r} not implemented in BioSTEAM (yet)')
-                    self.constraints = constraints = []
-                    stages = self.stages
-                    m, n = self.N_stages, self._N_chemicals
-                    last_stage = m - 1
-                    feed_flows, asplit_1, bsplit_1, _ = self._iter_args
-                    for i, stage in enumerate(stages):
-                        if i == 0:
-                            args = (i,)
-                            f = lambda x, i: feed_flows[i] - x[(i+1)*n:(i+2)*n] * asplit_1[i+1] - x[i*n:(i+1)*n] + 1e-6
-                        elif i == last_stage:
-                            args_last = args
-                            args = (i, f, args_last)
-                            f = lambda x, i, f, args_last: feed_flows[i] + f(x, *args_last) - x[i*n:] + 1e-6
-                        else:
-                            args_last = args
-                            args = (i, f, args_last)
-                            f = lambda x, i, f, args_last: feed_flows[i] + f(x, *args_last) - x[(i+1)*n:(i+2)*n] * asplit_1[i+1] - x[i*n:(i+1)*n] + 1e-6
-                        constraints.append(
-                            dict(type='ineq', fun=f, args=args)
-                        )
-                    result = solver(
-                        self._energy_balance_error_at_top_flow_rates, 
-                        self.get_top_flow_rates_flat(),
-                        constraints=constraints,
-                        bounds=[(0, None)] * (m * n),
-                        **options,
-                    )
-                    self.set_flow_rates(result.x.reshape([m, n]))
                 elif bounded and not constraints:
-                    # raise NotImplementedError(f'optimize method {self.method!r} not implemented in BioSTEAM (yet)')
                     self.iter = 0
                     partitions = self.partitions
                     Sb, safe = bottoms_stripping_factors_safe(
@@ -1754,8 +1694,8 @@ def _run(self):
                         bounds,
                         **options,
                     )
+                    if self._has_lle: self.update_energy_balance_temperatures()
                 elif not (constraints or bounded):
-                    # raise NotImplementedError(f'optimize method {self.method!r} not implemented in BioSTEAM (yet)')
                     self.iter = 0
                     partitions = self.partitions
                     Sb, safe = bottoms_stripping_factors_safe(
@@ -1785,11 +1725,12 @@ def _run(self):
                     self.set_flow_rates(
                         estimate_top_flow_rates(Sb, *self._iter_args, safe)
                     )
+                    if self._has_lle: self.update_energy_balance_temperatures()
             else:
                 raise RuntimeError(
                     f'invalid algorithm {algorithm!r}, only {self.available_algorithms} are allowed'
                 )
-            self.correct_overall_mass_balance()
+            # self.correct_overall_mass_balance()
         except Exception as e:
             if self.use_cache:
                 self.use_cache = False
@@ -2140,7 +2081,7 @@ def get_energy_balance_phase_ratio_departures(self):
                 hv[i] = top.h
             if Li == 0:
                 top.phase = 'l'
-                hl[i] = bottom.h
+                hl[i] = top.h
                 top.phase = 'g'
             else:
                 hl[i] = bottom.h
@@ -2171,7 +2112,7 @@ def get_energy_balance_phase_ratio_departures(self):
     def update_energy_balance_phase_ratio_departures(self):
         dBs = self.get_energy_balance_phase_ratio_departures()
         for i, dB in zip(self.partitions, dBs):
-            if i.B_specification is None: i.B += i.B_relaxation_factor * dB
+            if i.B_specification is None: i.B += (1 - i.B_relaxation_factor) * dB
     
     def update_energy_balance_temperatures(self):
         dTs = self.get_energy_balance_temperature_departures()
@@ -2182,6 +2123,7 @@ def update_energy_balance_temperatures(self):
             partition = stage.partition
             partition.T += dT
             for i in partition.outs: i.T += dT
+        return dTs
        
     def run_mass_balance(self):
         partitions = self.partitions
@@ -2291,17 +2233,34 @@ def _split_objective(self, splits):
         if not self._has_vle: raise NotImplementedError('only VLE objective function exists')
         self.set_flow_rates(top_flow_rates)
         P = self.P
-        for i in self.stages:
-            mixer = i.mixer
-            partition = i.partition
-            mixer.outs[0].mix_from(
-                mixer.ins, energy_balance=False,
+        stages = self.stages
+        if self._has_vle: 
+            for i in stages:
+                mixer = i.mixer
+                partition = i.partition
+                mixer.outs[0].mix_from(
+                    mixer.ins, energy_balance=False,
+                )
+                mixer.outs[0].P = P
+                partition._run_decoupled_KTvle(P=P, K_relaxation_factor=0, T_relaxation_factor=0)
+                partition._run_decoupled_B()
+                T = partition.T
+                for i in (partition.outs + i.outs): i.T = T
+            energy_error = lambda stage: abs(
+                (sum([i.H for i in stage.outs]) - sum([i.H for i in stage.ins])) / sum([i.C for i in stage.outs])
             )
-            mixer.outs[0].P = P
-            partition._run_decoupled_KTvle(P=P, K_relaxation_factor=0, T_relaxation_factor=0)
-            partition._run_decoupled_B()
-            T = partition.T
-            for i in (partition.outs + i.outs): i.T = T
+            total_energy_error = sum([energy_error(stage) for stage in stages if stage.B_specification is None and stage.T_specification is None])
+        else:
+            for i in range(5):
+                dTs = self.update_energy_balance_temperatures()
+                if np.abs(dTs).sum() < 1e-12: break
+            for i in self.stages: 
+                mixer = i.mixer
+                partition = i.partition
+                mixer.outs[0].mix_from(
+                    mixer.ins, energy_balance=False,
+                )
+                partition._run_lle(update=False, P=P)
         partitions = self.partitions
         if Sb_index:
             Sb_new = np.array([1 / (partitions[i].K * partitions[i].B) for i in Sb_index]).flatten()
@@ -2309,12 +2268,11 @@ def _split_objective(self, splits):
             Sb_new = np.array([1 / (i.K * i.B) for i in partitions]).flatten()
         splits_new = 1 / (Sb_new + 1)
         diff = splits_new - splits
-        energy_error = lambda stage: abs(
-            (sum([i.H for i in stage.outs]) - sum([i.H for i in stage.ins])) / sum([i.C for i in stage.outs])
-        )
-        total_energy_error = sum([energy_error(stage) for stage in self.stages if stage.B_specification is None and stage.T_specification is None])
         total_split_error = np.abs(diff).sum()
-        total_error = total_split_error + total_energy_error / self.N_stages
+        if self._has_vle:
+            total_error = total_split_error + total_energy_error / self.N_stages
+        else:
+            total_error = total_split_error
         err = np.sqrt(total_error)
         return err
     
@@ -2333,20 +2291,36 @@ def _lnSb_objective(self, lnSb):
             top_flow_rates = estimate_top_flow_rates(
                 Sb, *self._iter_args, True
             )
-        if not self._has_vle: raise NotImplementedError('only VLE objective function exists')
         self.set_flow_rates(top_flow_rates)
+        stages = self.stages
         P = self.P
-        for i in self.stages:
-            mixer = i.mixer
-            partition = i.partition
-            mixer.outs[0].mix_from(
-                mixer.ins, energy_balance=False,
+        if self._has_vle: 
+            for i in stages:
+                mixer = i.mixer
+                partition = i.partition
+                mixer.outs[0].mix_from(
+                    mixer.ins, energy_balance=False,
+                )
+                mixer.outs[0].P = P
+                partition._run_decoupled_KTvle(P=P, K_relaxation_factor=0, T_relaxation_factor=0)
+                partition._run_decoupled_B()
+                T = partition.T
+                for i in (partition.outs + i.outs): i.T = T
+            energy_error = lambda stage: abs(
+                (sum([i.H for i in stage.outs]) - sum([i.H for i in stage.ins])) / sum([i.C for i in stage.outs])
             )
-            mixer.outs[0].P = P
-            partition._run_decoupled_KTvle(P=P, K_relaxation_factor=0, T_relaxation_factor=0)
-            partition._run_decoupled_B()
-            T = partition.T
-            for i in (partition.outs + i.outs): i.T = T
+            total_energy_error = sum([energy_error(stage) for stage in stages if stage.B_specification is None and stage.T_specification is None])
+        else:
+            for i in range(5):
+                dTs = self.update_energy_balance_temperatures()
+                if np.abs(dTs).sum() < 1e-12: break
+            for i in self.stages: 
+                mixer = i.mixer
+                partition = i.partition
+                mixer.outs[0].mix_from(
+                    mixer.ins, energy_balance=False,
+                )
+                partition._run_lle(update=False, P=P)
         partitions = self.partitions
         if Sb_index:
             Sb_new = np.array([1 / (partitions[i].K * partitions[i].B) for i in Sb_index]).flatten()
@@ -2355,12 +2329,11 @@ def _lnSb_objective(self, lnSb):
         splits_new = 1 / (Sb_new + 1)
         splits = 1 / (Sb_original + 1)
         diff = splits_new - splits
-        energy_error = lambda stage: abs(
-            (sum([i.H for i in stage.outs]) - sum([i.H for i in stage.ins])) / sum([i.C for i in stage.outs])
-        )
-        total_energy_error = sum([energy_error(stage) for stage in self.stages if stage.B_specification is None and stage.T_specification is None])
         total_split_error = np.abs(diff).sum()
-        total_error = total_split_error + total_energy_error / self.N_stages
+        if self._has_vle:
+            total_error = total_split_error + total_energy_error / self.N_stages
+        else:
+            total_error = total_split_error
         err = np.sqrt(total_error)
         return err
     
@@ -2382,7 +2355,8 @@ def _iter(self, top_flow_rates):
                     partition._run_decoupled_KTvle(P=P)
                     T = partition.T
                     for i in (partition.outs + i.outs): i.T = T
-                    if partition.reaction: partition._run_decoupled_reaction(P=P)
+                    if partition.reaction: 
+                        partition._run_decoupled_reaction(P=P)
             else:
                 for i in stages:
                     mixer = i.mixer
@@ -2411,11 +2385,9 @@ def psuedo_equilibrium(top_flow_rates):
                             partition._run_decoupled_Kgamma(P=P)
                         self.interpolate_missing_variables()
                         return self.run_mass_balance()
-                options = tmo.equilibrium.LLE.pseudo_equilibrium_inner_loop_options.copy()
-                options['xtol'] *= self.feed_flows.max()
                 self.set_flow_rates(
                     flx.fixed_point(
-                        psuedo_equilibrium, top_flow_rates, **options,
+                        psuedo_equilibrium, top_flow_rates, **self._psuedo_equilibrium_options,
                     )
                 )
             for i in stages: 
@@ -2531,7 +2503,7 @@ def _conditional_iter(self, top_flow_rates):
         return top_flow_rates_new, not_converged
 
     def _sequential_iter(self, top_flow_rates):
-        self.fallback_iter += 1
+        self.iter += 1
         self.set_flow_rates(top_flow_rates)
         for i in self.stages: i._run()
         for i in reversed(self.stages): i._run()
@@ -2547,7 +2519,7 @@ def _sequential_iter(self, top_flow_rates):
                 max_errors = np.maximum.reduce([abs(mol[nonzero_index]), abs(mol_new[nonzero_index])])
                 rmol_error = (mol_errors / max_errors).max()
                 not_converged = (
-                    self.fallback_iter < self.fallback_maxiter and (mol_error > self.molar_tolerance
+                    self.iter < self.maxiter and (mol_error > self.molar_tolerance
                      or rmol_error > self.relative_molar_tolerance)
                 )
             else:
diff --git a/tests/test_equilibrium_stages.py b/tests/test_equilibrium_stages.py
index 122b887f..4c681f3f 100644
--- a/tests/test_equilibrium_stages.py
+++ b/tests/test_equilibrium_stages.py
@@ -3,6 +3,10 @@
 """
 import pytest
 import biosteam as bst
+import thermosteam as tmo
+from math import exp, log
+import numpy as np
+from thermosteam.constants import R
 from numpy.testing import assert_allclose
 
 def test_multi_stage_adiabatic_vle():
@@ -117,7 +121,7 @@ def test_multi_stage_adiabatic_vle():
         rtol=0.01,
     )
    
-def test_reactive_distillation():
+def test_lactic_acid_ethanol_reactive_distillation():
     import biosteam as bst
     from thermosteam.constants import R
     from math import exp
@@ -137,7 +141,7 @@ def rate(self, stream):
             LaEt, La, H2O, EtOH = stream.mol / stream.F_mol
             return 3600 * (kf * La * EtOH - kr * LaEt * H2O) # kmol / kg-catalyst / hr
     
-    rxn = Esterification('LacticAcid + Ethanol -> H2O + EthylLactate')
+    rxn = Esterification('LacticAcid + Ethanol -> H2O + EthylLactate', reactant='LacticAcid')
     stream = bst.Stream(
         H2O=10, Ethanol=10, LacticAcid=2, T=355,
     )
@@ -158,10 +162,155 @@ def rate(self, stream):
     )
     distillation.simulate()
     flows = [
-        [6.575895653800856e-06, 3.6651391022703107e-09, 
-         3.1210737204600005, 8.043951319613521],
-        [0.019035234263940305, 2.0190418064944553, 
-         6.897968089699594, 1.9750904905460722],
+        [6.585111268569887e-06, 3.6465771244802687e-09,
+         3.110437134095077, 8.015523262914513],
+        [0.019007540635984595, 1.9809858706061696,
+         6.908576991652181, 1.9654626113382356],
+    ]
+    for i, j in zip(distillation.outs, flows):    
+        assert_allclose(i.mol, j, rtol=1e-5, atol=1e-3)
+    sys = bst.System.from_units(units=[distillation])
+    sys._setup()
+    sys.run_phenomena()
+    for i, j in zip(distillation.outs, flows):    
+        assert_allclose(i.mol, j, rtol=1e-5, atol=1e-3)
+    
+def test_acetic_acid_reactive_distillation():
+    import biosteam as bst
+    import thermosteam as tmo
+    from thermosteam.constants import R
+    import numpy as np
+    from math import exp
+    from numpy.testing import assert_allclose
+    bst.settings.set_thermo(['Water', 'AceticAcid', 'MethylAcetate', 'Methanol'], cache=True)
+    Gamma = tmo.equilibrium.DortmundActivityCoefficients(bst.settings.chemicals)
+    class Esterification(bst.KineticReaction):
+        def volume(self, stream): # kg of catalyst
+            rho_cat = 770 # kg / m3
+            liquid_volume = self.liquid_volume
+            catalyst_volume = 0.5 * liquid_volume
+            catalyst_mass = catalyst_volume * rho_cat
+            return catalyst_mass
+        
+        def rate(self, stream): # kmol/kg-catalyst/hr
+            F_mol = stream.F_mol
+            if not F_mol: return 0
+            K_AcOH = 3.15
+            K_MeOH = 5.64
+            K_MeOAc = 4.15
+            K_H2O = 5.24
+            K_adsorption = np.array([K_H2O, K_AcOH, K_MeOAc, K_MeOH])
+            MWs = stream.chemicals.MW
+            k0f = 8.497e6
+            EAf = 60.47
+            k0r = 6.127e5
+            EAr = 63.73
+            T = stream.T
+            mol = stream.mol
+            x = mol / F_mol
+            gamma = Gamma(x, T)
+            activity = gamma * x
+            a = K_adsorption * activity / MWs
+            H2O, LA, BuLA, BuOH = activity
+            kf = k0f * exp(-EAf / (R * T))
+            kr = k0r * exp(-EAr / (R * T))
+            a_sum = sum(a)
+            return 3600 * (kf * LA * BuOH - kr * BuLA * H2O) / (a_sum * a_sum) # kmol / kg-catalyst / hr
+        
+    methanol = bst.Stream(
+        Methanol=30,
+    )
+    methanol.vle(V=1, P=101325)
+    acetic_acid = bst.Stream(
+        AceticAcid=10, Water=200
+    )
+    acetic_acid.vle(V=0, P=101325)
+    stage_reactions = {
+        i: Esterification('AceticAcid + Methanol -> MethylAcetate + Water', reactant='AceticAcid')
+        for i in range(3, 9)
+    }
+    distillation = bst.MESHDistillation(
+        N_stages=10,
+        ins=[methanol, acetic_acid],
+        feed_stages=[3, 8],
+        outs=['distillate', 'bottoms_product'],
+        full_condenser=False,
+        reflux=1.0,
+        boilup=2.0,
+        use_cache=True,
+        LHK=('MethylAcetate', 'Water'),
+        stage_reactions=stage_reactions,
+        method='fixed-point',
+        maxiter=100,
+    )
+    distillation.simulate()
+    flows = [
+        [6.585111268569887e-06, 3.6465771244802687e-09,
+         3.110437134095077, 8.015523262914513],
+        [0.019007540635984595, 1.9809858706061696,
+         6.908576991652181, 1.9654626113382356],
+    ]
+    for i, j in zip(distillation.outs, flows):    
+        assert_allclose(i.mol, j, rtol=1e-5, atol=1e-3)
+    sys = bst.System.from_units(units=[distillation])
+    sys._setup()
+    sys.run_phenomena()
+    for i, j in zip(distillation.outs, flows):    
+        assert_allclose(i.mol, j, rtol=1e-5, atol=1e-3)
+    
+    
+def test_lactic_acid_butanol_reactive_distillation():
+    import biosteam as bst
+    from thermosteam.constants import R
+    from math import exp
+    from numpy.testing import assert_allclose
+    bst.settings.set_thermo(['Water', 'LacticAcid', 'ButylLactate', 'Butanol'], cache=True)
+    
+    class Esterification(bst.KineticReaction):
+        
+        def volume(self, stream): # kg of catalyst
+            rho_cat = 770 # kg / m3
+            liquid_volume = self.liquid_volume
+            catalyst_volume = 0.5 * liquid_volume
+            catalyst_mass = catalyst_volume * rho_cat
+            return catalyst_mass
+        
+        def rate(self, stream): # kmol/kg-catalyst/hr
+            T = stream.T
+            # if T > 370: return 0 # Prevents multiple steady states
+            kf = 2.59e4 * exp(-5.340e4 / (R * T))
+            kr = 3.80e3 * exp(-5.224e4 / (R * T))
+            H2O, LA, BuLA, BuOH = stream.mol / stream.F_mol
+            return 3600 * (kf * LA * BuOH - kr * BuLA * H2O) # kmol / kg-catalyst / hr
+        
+    stream = bst.Stream(
+        H2O=10, Butanol=10, LacticAcid=2, T=355,
+    )
+    stream.vle(V=0, P=101325)
+    stage_reactions = {
+        i: Esterification('LacticAcid + Butanol -> Water + ButylLactate', reactant='LacticAcid')
+        for i in range(2)
+    }
+    distillation = bst.MESHDistillation(
+        N_stages=5,
+        ins=[stream],
+        feed_stages=[2],
+        outs=['distillate', 'bottoms_product'],
+        full_condenser=False,
+        reflux=1.0,
+        boilup=2.0,
+        use_cache=True,
+        LHK=('Butanol', 'ButylLactate'),
+        stage_reactions=stage_reactions,
+        method='fixed-point',
+        maxiter=100,
+    )
+    distillation.simulate()
+    flows = [
+        [6.585111268569887e-06, 3.6465771244802687e-09,
+         3.110437134095077, 8.015523262914513],
+        [0.019007540635984595, 1.9809858706061696,
+         6.908576991652181, 1.9654626113382356],
     ]
     for i, j in zip(distillation.outs, flows):    
         assert_allclose(i.mol, j, rtol=1e-5, atol=1e-3)
@@ -211,7 +360,7 @@ def test_distillation():
 if __name__ == '__main__':
     test_multi_stage_adiabatic_vle()
     test_distillation()
-    test_reactive_distillation()
+    test_lactic_acid_ethanol_reactive_distillation()
 
 
 
diff --git a/thermosteam b/thermosteam
index 800b77fb..b6e4e292 160000
--- a/thermosteam
+++ b/thermosteam
@@ -1 +1 @@
-Subproject commit 800b77fb823198274f179b95deec9bc11683b915
+Subproject commit b6e4e292568787cf3cb41f3e0c2290592cb2336b