Skip to content

Commit

Permalink
better PO-alg
Browse files Browse the repository at this point in the history
  • Loading branch information
yoelcortes committed Aug 2, 2024
1 parent 907714c commit 8aeee1c
Show file tree
Hide file tree
Showing 8 changed files with 389 additions and 223 deletions.
2 changes: 1 addition & 1 deletion benchmark/butanol_purification_system.py
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,7 @@ def register(*args, **kwargs):

@register(
'butanol_purification', 'Butanol purification',
4, [0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4], 'BtOH\nsep.'
0.5, [0.1, 0.2, 0.3, 0.4, 0.5], 'BtOH\nsep.'
)
def create_system_butanol_purification(alg):
bst.settings.set_thermo(['Water', 'Butanol'], cache=True)
Expand Down
2 changes: 1 addition & 1 deletion benchmark/ethanol_purification_system.py
Original file line number Diff line number Diff line change
Expand Up @@ -13,7 +13,7 @@

@register(
'ethanol_purification', 'Ethanol purification',
0.1, [0.02, 0.04, 0.06, 0.08, 0.01], 'EtOH\nsep.',
0.15, [0.03, 0.06, 0.09, 0.12, 0.15], 'EtOH\nsep.',
)
def create_system_ethanol_purification(alg):
bst.settings.set_thermo(['Water', 'Ethanol'], cache=True)
Expand Down
149 changes: 149 additions & 0 deletions benchmark/haber_bosch_process.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,149 @@
# -*- 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_haber_bosch_process',
)

@register(
'haber_bosch_process', 'Haber-Bosch process',
10, [2, 4, 6, 8, 10], 'HB\nprocess.'
)
def create_system_haber_bosch_process(alg='sequential modular'):
bst.settings.set_thermo(['N2', 'H2', 'NH3'], cache=True)
bst.settings.mixture.include_excess_energies = True

with bst.System(algorithm=alg) as sys:
feed = bst.Stream(
'feed',
H2=3,
N2=1,
P=300e5,
phase='g'
)
preheater = bst.HXutility(ins=feed, T=450 + 273.15)
reactor = bst.
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()
27 changes: 12 additions & 15 deletions biosteam/_system.py
Original file line number Diff line number Diff line change
Expand Up @@ -142,8 +142,6 @@ def _solve_material_flows(self, composition_sensitive):
A = []
b = []
for node in nodes:
if not node._create_material_balance_equations:
raise NotImplementedError(f'{node!r} has no method `_create_material_balance_equations`')
for coefficients, value in node._create_material_balance_equations(composition_sensitive):
coefficients = {
streams[i.imol]: j for i, j in coefficients.items()
Expand Down Expand Up @@ -190,17 +188,19 @@ def solve_energy_departures(self):
A = []
b = []
for node in nodes:
if not node._create_energy_departure_equations:
raise NotImplementedError(f'{node!r} has no method `_create_energy_departure_equations`')
for coefficients, value in node._create_energy_departure_equations():
A.append(coefficients)
b.append(value)
A, objs = dictionaries2array(A)
departures = solve(A, np.array(b).T).T
for obj, departure in zip(objs, departures):
if not obj._update_energy_variable:
try:
for obj, departure in zip(objs, departures):
obj._update_energy_variable(departure)
except AttributeError as e:
if obj._update_energy_variable:
raise e
else:
raise NotImplementedError(f'{obj!r} has no method `_update_energy_variable`')
obj._update_energy_variable(departure)

def __enter__(self):
units = self.stages
Expand Down Expand Up @@ -2363,7 +2363,6 @@ def stage_configuration(self, aggregated=False):
else:
return self._stage_configuration
except:
feeds = [i for i in self.feeds if i.material_equations]
stages = self.stages
streams = list(set([i for s in stages for i in s.ins + s.outs]))
connections = [(i, i.source, i.sink) for i in streams]
Expand All @@ -2387,21 +2386,19 @@ def run_phenomena(self):
"""Decouple and linearize material, equilibrium, summation, enthalpy,
and reaction phenomena and iteratively solve them."""
path = self.unit_path
n = -1
with self.stage_configuration(aggregated=False) as conf:
try:
conf.solve_nonlinearities()
conf.solve_energy_departures()
conf.solve_material_flows()
for n, i in enumerate(path):
for i in path:
i.run()
conf.solve_energy_departures()
conf.solve_material_flows()
except Exception as e:
if self._iter > 0:
raise e
breakpoint()
for i in path[n+1:]: i.run()
except NotImplementedError as error:
raise error
except:
for i in path: i.run()

# try:
# with self.stage_configuration(aggregated=True) as conf:
Expand Down
44 changes: 31 additions & 13 deletions biosteam/_unit.py
Original file line number Diff line number Diff line change
Expand Up @@ -253,29 +253,47 @@ def __init_subclass__(cls,

def _begin_equations(self, composition_sensitive):
inlets = self.ins
if composition_sensitive:
material_equations = [f() for i in inlets for f in i.material_equations]
fresh_inlets = []
process_inlets = []
material_equations = [f() for i in inlets for f in i.material_equations]
if not material_equations:
if composition_sensitive:
for i in inlets:
if i.isfeed() or not getattr(i.source, 'composition_sensitive', False):
fresh_inlets.append(i)
else:
process_inlets.append(i)
else:
for i in inlets:
if i.isfeed():
fresh_inlets.append(i)
else:
process_inlets.append(i)
equations = material_equations
elif composition_sensitive:
equations = []
dependent_streams = []
independent_streams = []
for eq in material_equations:
dct, value = eq
for i in dct:
exclude = i.source and not getattr(i.source, 'composition_sensitive', False)
if exclude:
independent_streams.append(i.imol)
break
if i.source and not getattr(i.source, 'composition_sensitive', False): break
else:
dependent_streams.append(i.imol)
equations.append(eq)
independent_streams = set(independent_streams)
dependent_streams = set(dependent_streams)
fresh_inlets = [i for i in inlets if i.imol in independent_streams or i.isfeed() and not i.imol in dependent_streams]
process_inlets = [i for i in inlets if (not i.isfeed() and i.imol not in independent_streams) or i.imol in dependent_streams]
for i in inlets:
isfeed = i.isfeed() or not getattr(i.source, 'composition_sensitive', False)
if i.imol not in dependent_streams and isfeed:
fresh_inlets.append(i)
else:
process_inlets.append(i)
else:
fresh_inlets = [i for i in inlets if i.isfeed() and not i.material_equations]
process_inlets = [i for i in inlets if not i.isfeed() or i.material_equations]
equations = [f() for i in inlets for f in i.material_equations]
for i in inlets:
if i.isfeed() and not i.material_equations:
fresh_inlets.append(i)
else:
process_inlets.append(i)
equations = material_equations
return fresh_inlets, process_inlets, equations

Inlets = piping.Inlets
Expand Down
8 changes: 5 additions & 3 deletions biosteam/units/stage.py
Original file line number Diff line number Diff line change
Expand Up @@ -602,7 +602,8 @@ def _update_equilibrium_variables(self):
T = partition.T
for i in (partition.outs, self.outs): i.T = T
elif phases == ('L', 'l'):
self.partition._run_lle(single_loop=True)
# self.partition._run_lle(single_loop=True)
pass
else:
raise NotImplementedError(f'K for phases {phases} is not yet implemented')

Expand Down Expand Up @@ -649,7 +650,7 @@ class PhasePartition(Unit):
strict_infeasibility_check = False
dmol_relaxation_factor = 0.9
S_relaxation_factor = 0
B_relaxation_factor = 0.6
B_relaxation_factor = 0.5

def _init(self, phases, partition_data, top_chemical=None, reaction=None):
self.partition_data = partition_data
Expand Down Expand Up @@ -1528,7 +1529,8 @@ def _update_nonlinearities(self):
if self._has_vle:
for i in self.stages: i._update_nonlinearities()
elif self._has_lle:
self.update_lle_variables()
pass
# self.update_lle_variables()

def _update_aggretated_nonlinearities(self):
top_flow_rates = self.get_top_flow_rates()
Expand Down
Loading

0 comments on commit 8aeee1c

Please sign in to comment.