even more robust PO-alg

biosteam/_system.py

@@ -79,7 +79,9 @@ def __call__(self): self.f(*self.args)
 class Configuration:
  __slots__ = ('path', 'stages', 'nodes', 'streams',
  'stream_ref', 'connections', 'aggregated',
- 'composition_sensitive_path', '_has_dynamic_coefficients')
+ 'composition_sensitive_path', 
+ 'composition_sensitive_nodes',
+ '_has_dynamic_coefficients')
 
  def __init__(self, path, stages, nodes, streams, stream_ref, connections, aggregated):
  self.path = path
@@ -90,6 +92,7 @@ def __init__(self, path, stages, nodes, streams, stream_ref, connections, aggreg
  self.connections = connections
  self.aggregated = aggregated
  self.composition_sensitive_path = [i for i in path if getattr(i, 'composition_sensitive', False)]
+ self.composition_sensitive_nodes = [i for i in nodes if getattr(i, 'composition_sensitive', False) or isinstance(i, Stream)]
 
  def solve_nonlinearities(self):
  if self.aggregated:
@@ -106,20 +109,19 @@ def solve_nonlinearities(self):
  def solve_material_flows(self):
  if self.composition_sensitive_path:
  for i in self.composition_sensitive_path: i._update_composition_parameters()
- flows = self._solve_material_flows()
+ flows = self._solve_material_flows(composition_sensitive=True)
 
  def update_inner_material_balance_parameters(flows):
- self.solve_nonlinearities
  for i in self.composition_sensitive_path: i._update_composition_parameters()
- return self._solve_material_flows()
+ return self._solve_material_flows(composition_sensitive=True)
 
  flx.fixed_point(
  update_inner_material_balance_parameters,
  flows, xtol=1e-6 + 1e-6 * flows.max(), maxiter=3, checkiter=False,
  checkconvergence=False,
  )
  for i in self.composition_sensitive_path: i._update_net_flow_parameters()
- self._solve_material_flows()
+ self._solve_material_flows(composition_sensitive=False)
 
  def dynamic_coefficients(self, b):
  try:
@@ -131,15 +133,18 @@ def dynamic_coefficients(self, b):
  delayed = [(i, j) for i, j in enumerate(b) if j.dtype is ObjectType]
  return delayed
 
- def _solve_material_flows(self):
- nodes = self.nodes
+ def _solve_material_flows(self, composition_sensitive):
+ if composition_sensitive:
+ nodes = self.composition_sensitive_nodes
+ else:
+ nodes = self.nodes
  streams = self.stream_ref
  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():
+ for coefficients, value in node._create_material_balance_equations(composition_sensitive):
  coefficients = {
  streams[i.imol]: j for i, j in coefficients.items()
  }
@@ -2364,17 +2369,15 @@ def stage_configuration(self, aggregated=False):
  connections = [(i, i.source, i.sink) for i in streams]
  if aggregated:
  stages = self.aggregated_stages
- nodes = stages + feeds
  streams = [i for u in stages for i in u.ins + u.outs]
  stream_ref = {i.imol: i for u in stages for i in (*u.ins, *u.outs)}
- nodes = [i for i in nodes if not getattr(i, 'decoupled', False)]
+ nodes = [i for i in stages if not getattr(i, 'decoupled', False)]
  self._aggregated_stage_configuration = conf = Configuration(
  self.path, stages, nodes, streams, stream_ref, connections, aggregated
  )
  else:
- nodes = stages + feeds
  stream_ref = {i.imol: i for i in streams}
- nodes = [i for i in nodes if not getattr(i, 'decoupled', False)]
+ nodes = [i for i in stages if not getattr(i, 'decoupled', False)]
  self._stage_configuration = conf = Configuration(
  self.path, stages, nodes, streams, stream_ref, connections, aggregated
  )
@@ -2387,16 +2390,19 @@ def run_phenomena(self):
  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):
  i.run()
  conf.solve_energy_departures()
  conf.solve_material_flows()
- except: 
+ except Exception as e: 
+ if self._iter > 0: 
+ raise e
+ breakpoint() 
  for i in path[n+1:]: i.run()
- else:
- conf.solve_nonlinearities()
- conf.solve_energy_departures()
- conf.solve_material_flows()
+
  # try:
  # with self.stage_configuration(aggregated=True) as conf:
  # conf.solve_nonlinearities()

biosteam/_unit.py

@@ -251,6 +251,33 @@ def __init_subclass__(cls,
  #: Update reaction conversion for phenomena-oriented simulation.
  _update_reaction_conversion = AbstractMethod
 
+ 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]
+ 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
+ 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]
+ 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]
+ return fresh_inlets, process_inlets, equations
+
  Inlets = piping.Inlets
  Outlets = piping.Outlets
 

biosteam/units/compressor.py

@@ -568,7 +568,8 @@ def _design(self):
  super()._design()
  self._set_power(self.design_results['Ideal power'] / self.eta)
 
- def _create_energy_departure_equations(self):
+ def _create_material_balance_equations(self, composition_sensitive):
+ fresh_inlets, process_inlets, equations = self._begin_equations(composition_sensitive)
  # Ll: C1dT1 - Ce2*dT2 - Cr0*dT0 - hv2*L2*dB2 = Q1 - H_out + H_in
  # gl: hV1*L1*dB1 - hv2*L2*dB2 - Ce2*dT2 - Cr0*dT0 = Q1 + H_in - H_out
  feed = self.ins[0]
@@ -579,25 +580,7 @@ def _create_energy_departure_equations(self):
  J = R / (2. * product.Cn * self.eta) * log(feed.P / self.P)
  coeff = {self: (1 - J) / (1 + J) * product.Cn}
  for i in self.ins: i._update_energy_departure_coefficient(coeff)
- return [(coeff, 0)]
-
- def _create_material_balance_equations(self):
- inlets = self.ins
- fresh_inlets = [i for i in inlets if i.isfeed() and not i.equations]
- process_inlets = [i for i in inlets if not i.isfeed() or i.equations]
- product, = self.outs
- ones = np.ones(self.chemicals.size)
- minus_ones = -ones
- zeros = np.zeros(self.chemicals.size)
-
- # Overall flows
- eq_overall = {}
- for i in self.outs: eq_overall[i] = ones
- for i in process_inlets: eq_overall[i] = minus_ones
- equations = [
- (eq_overall, sum([i.mol for i in fresh_inlets], zeros))
- ]
-
+ equations.append((coeff, 0))
  return equations
 
  def _update_energy_variable(self, departure):

biosteam/units/mixing.py

@@ -123,11 +123,9 @@ def _create_energy_departure_equations(self):
  for i in self.ins: i._update_energy_departure_coefficient(coeff)
  return [(coeff, self.H_in - self.H_out)]
 
- def _create_material_balance_equations(self):
- inlets = self.ins
+ def _create_material_balance_equations(self, composition_sensitive):
+ fresh_inlets, process_inlets, equations = self._begin_equations(composition_sensitive)
  outlet, = self.outs
- 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]
  if len(outlet) == 1:
  ones = np.ones(self.chemicals.size)
  minus_ones = -ones
@@ -136,12 +134,11 @@ def _create_material_balance_equations(self):
  # Overall flows
  eq_overall = {outlet: ones}
  for i in process_inlets: eq_overall[i] = minus_ones
- equations = [
+ equations.append(
  (eq_overall, sum([i.mol for i in fresh_inlets], zeros))
- ]
+ )
  else:
  top, bottom = outlet
- equations = []
  ones = np.ones(self.chemicals.size)
  minus_ones = -ones
  zeros = np.zeros(self.chemicals.size)

biosteam/units/splitting.py

@@ -301,15 +301,12 @@ def _run(self):
  if self.P:
  top.P = bottom.P = self.P
 
- def _create_material_balance_equations(self):
- inlets = self.ins
+ def _create_material_balance_equations(self, composition_sensitive):
+ fresh_inlets, process_inlets, equations = self._begin_equations(composition_sensitive)
  top, bottom = self.outs
- equations = []
  ones = np.ones(self.chemicals.size)
  minus_ones = -ones
  zeros = np.zeros(self.chemicals.size)
- 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]
 
  # Overall flows
  eq_overall = {}

biosteam/units/stage.py

@@ -110,11 +110,9 @@ def _create_energy_departure_equations(self):
  for i in self.ins: i._update_energy_departure_coefficient(coeff)
  return [(coeff, outlet.H - sum([i.H for i in self.ins]))]
 
- def _create_material_balance_equations(self):
+ def _create_material_balance_equations(self, composition_sensitive):
  outlet = self.outs[0]
- inlets = self.ins
- 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]
+ fresh_inlets, process_inlets, equations = self._begin_equations(composition_sensitive)
  ones = np.ones(self.chemicals.size)
  minus_ones = -ones
  zeros = np.zeros(self.chemicals.size)
@@ -126,9 +124,10 @@ def _create_material_balance_equations(self):
  del eq_overall[i]
  else:
  eq_overall[i] = minus_ones
- return [
+ equations.append(
  (eq_overall, sum([i.mol for i in fresh_inlets], zeros))
- ]
+ )
+ return equations
 
  def _update_energy_variable(self, value):
  if self.T is None: return
@@ -160,15 +159,12 @@ def _run(self):
  self.reaction(outlet, Q=self.Q)
  outlet.T = self.T
 
- def _create_material_balance_equations(self):
- inlets = self.ins
+ def _create_material_balance_equations(self, composition_sensitive=False):
  product, = self.outs
- equations = []
  n = self.chemicals.size
  ones = np.ones(n)
  minus_ones = -ones
- 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]
+ fresh_inlets, process_inlets, equations = self._begin_equations(composition_sensitive)
  reaction = self.reaction
  if isinstance(reaction, bst.KRxn):
  # Overall flows
@@ -236,8 +232,6 @@ class StageEquilibrium(Unit):
  _ins_size_is_fixed = False
  _outs_size_is_fixed = False
  auxiliary_unit_names = ('partition', 'mixer', 'splitters')
- S_relaxation_factor = 0
- B_relaxation_factor = 0.5
 
  def __init__(self, ID='', ins=None, outs=(), thermo=None, *, 
  phases, partition_data=None, top_split=0, bottom_split=0,
@@ -460,7 +454,7 @@ def _create_energy_departure_equations(self):
  dH = self.Q + self.H_in - self.H_out
  return [(coeff, dH)]
 
- def _create_material_balance_equations(self):
+ def _create_material_balance_equations(self, composition_sensitive):
  phases = self.phases
  partition = self.partition
  chemicals = self.chemicals
@@ -484,13 +478,10 @@ def _create_material_balance_equations(self):
  self._update_separation_factors()
  top_split = self.top_split
  bottom_split = self.bottom_split
- inlets = self.ins
- 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]
+ fresh_inlets, process_inlets, equations = self._begin_equations(composition_sensitive)
  top, bottom, *_ = self.outs
  top_side_draw = self.top_side_draw
  bottom_side_draw = self.bottom_side_draw
- equations = []
  N = self.chemicals.size
  ones = np.ones(N)
  minus_ones = -ones
@@ -586,7 +577,7 @@ def _update_energy_variable(self, departure):
  top, bottom = partition.outs
  IDs = partition.IDs
  B = top.imol[IDs].sum() / bottom.imol[IDs].sum()
- self.B = B + (1 - self.B_relaxation_factor) * departure
+ self.B = 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
@@ -646,8 +637,8 @@ def _update_separation_factors(self, f=None):
  if self.phases == ('L', 'l'):
  self.Sb = Sb
  else:
- if f is None: f = self.S_relaxation_factor
- self.Sb = (1 - f) * Sb + f * self.Sb
+ if f is None: f = self.partition.S_relaxation_factor
+ self.Sb = (1 - f) * Sb + f * self.Sb if f else Sb
 
 
 # %%
@@ -657,6 +648,8 @@ class PhasePartition(Unit):
  _N_outs = 2
  strict_infeasibility_check = False
  dmol_relaxation_factor = 0.9
+ S_relaxation_factor = 0
+ B_relaxation_factor = 0.6
 
  def _init(self, phases, partition_data, top_chemical=None, reaction=None):
  self.partition_data = partition_data
@@ -1101,7 +1094,7 @@ class MultiStageEquilibrium(Unit):
  _N_ins = 2
  _N_outs = 2
  max_attempts = 10
- default_maxiter = 30
+ default_maxiter = 20
  default_fallback = None
  default_molar_tolerance = 1e-3
  default_relative_molar_tolerance = 1e-6
@@ -1460,7 +1453,7 @@ def _create_energy_departure_equations(self):
  else:
  return [(coeff, self.H_in - self.H_out + sum([i.Q for i in self.stages]))]
 
- def _create_material_balance_equations(self):
+ def _create_material_balance_equations(self, composition_sensitive):
  top, bottom = self.outs
  try:
  B = self.B
@@ -1483,11 +1476,8 @@ def _create_material_balance_equations(self):
  self.K = K = y / x
  self.B = B = F_top / F_bottom
 
- inlets = self.ins
- 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]
+ fresh_inlets, process_inlets, equations = self._begin_equations(composition_sensitive)
  top, bottom, *_ = self.outs
- equations = []
  ones = np.ones(self.chemicals.size)
  minus_ones = -ones
  zeros = np.zeros(self.chemicals.size)
@@ -1700,16 +1690,19 @@ def _run(self):
  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-1):
+ for n in range(self.max_attempts-1):
  self.iter = 0
- self.attempt = i
+ self.attempt = n
  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()
  top_flow_rates = self.get_top_flow_rates()
+ # for i in self.partitions: i.B_relaxation_factor += (0.9 - i.B_relaxation_factor) * 0.5
  else:
  break
+ # if n:
+ # for i in self.partitions: del i.B_relaxation_factor
  if self.iter == self.maxiter:
  top_flow_rates = solver(self._conditional_iter, top_flow_rates)
  if self.iter == self.maxiter and self.fallback and self.fallback != (self.algorithm, self.method):
@@ -2173,9 +2166,8 @@ def get_energy_balance_phase_ratio_departures(self):
 
  def update_energy_balance_phase_ratio_departures(self):
  dBs = self.get_energy_balance_phase_ratio_departures()
- g = (1 - StageEquilibrium.B_relaxation_factor)
  for i, dB in zip(self.partitions, dBs):
- if i.B_specification is None: i.B += g * 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()