Optimizer Implementation in SMB

bpimentel · April 8, 2024, 5:40am

New user here - really struggling to understand how to execute an optimization problem in CADET

I’m using the Carousel Builder example as is from the documentation as a starting point, but I wanted to add an optimizer for, lets say, eluent minimization.

I did add the stationary evaluator to the program

from CADETProcess.simulator import Cadet
from CADETProcess.stationarity import StationarityEvaluator
from CADETProcess.stationarity import RelativeArea

process_simulator = Cadet()

evaluator = StationarityEvaluator()

criterion = RelativeArea()
criterion.threshold = 5e-3
evaluator.add_criterion(criterion)

process_simulator.stationarity_evaluator = evaluator
process_simulator.evaluate_stationarity = True

This is what I have so far - allowing extract and raffinate draws to vary, along with eluent flow rate and switch time, keeping feed inlet constant.

I know I would still need to add a purity and recovery constraint to make it all work.

optimization_problem = OptimizationProblem('eluent_usage')
optimization_problem.add_variable('w_e', lb=0, ub=1)
optimization_problem.add_variable('w_r', lb=0, ub=1)
optimization_problem.add_variable('eluent.flow_rate', lb=1e-7, ub=10e-7)
optimization_problem.add_variable('builder.switch_time', lb = 150, ub = 600)

The general optimization cases I’d like to run would be

Eluent minimization : obj: Eluent/Feed ----- at a minimum purity and recovery, constant feed
Productivity: obj: Component flow ----- at raffinate/extract at a minimum purity and recovery, constant feed
Recovery: obj: Component out/Component in ------ at constant feed, minimum purity
Purity: obj: Average component purity at raff/extract — at constant feed, minimum recovery.
This last one I at least know how to access the variable after the fact, but don’t know if its accessible during the simulation.

simulation_results.solution_cycles.extract.inlet[-1].local_purity_species.mean(axis=0)

Are the chromatography functions like “Productivity” usable in a carousel mode?

Additionally - any idea how long these optimizations tend to take? I think i hit CSS in 120 s.

j.schmoelder · April 8, 2024, 4:24pm

Hi Brian and welcome to the CADET-Forum!

Optimizing SMB systems is certainly not trivial but it can be done. It involves some advanced features of CADET-Process which might seem a bit unintuitive at first but I haven’t found a more elegant solution yet. So I’m interested in your feedback.

First of all, calculating key performance indicators such as eluent consumption, recovery, purity etc. requires some additional post-processing in CADET-Process. For this purpose, we provide a fractionation module.

Depending on how you want to set up your optimization problem, this can be trivial, because you simply define a single fraction per outlet, that simply collects the entire stream. However, since we also allow automatic determination of fractionation times, I would suggest setting up a FractionationOptimizer.

You can include this in your processing pipeline for defining optimization problems. For a very similar example, refer to one of our SMB case studies.

process = smb_builder.build_process()

optimization_problem.add_evaluation_object(process)

optimization_problem.add_variable(
    'switch_time.time',
    lb=150, ub=600,
    transform='auto'
)
optimization_problem.add_variable(
    'cycle_time',
)
optimization_problem.add_variable_dependency(
    'cycle_time', 'switch_time.time', lambda x: 4 * x
)

optimization_problem.add_variable(
    'flow_sheet.eluent.flow_rate',
    lb=1e-7, ub=10e-7,
    transform='auto'
  )

optimization_problem.add_variable(
    name='w_r',
    parameter_path='flow_sheet.output_states.zone_III_outlet',
    lb=0, ub=1,
    pre_processing=lambda x: [x[0], 1 - x[0]],
)
optimization_problem.add_variable(
    name='w_e',
    parameter_path='flow_sheet.output_states.zone_I_outlet',
    lb=0, ub=1,
    pre_processing=lambda x: [x[0], 1 - x[0]],
)

# Setup Simulator
process_simulator = Cadet(options.cadet_path)
process_simulator.evaluate_stationarity = True

optimization_problem.add_evaluator(process_simulator)

# Setup Fractionator
frac_opt = FractionationOptimizer()

optimization_problem.add_evaluator(
    frac_opt,
    kwargs={'purity_required': [0.95, 0.95]}
)

# Setup Objectives
eluent_consumption = EluentConsumption()
optimization_problem.add_objective(
    eluent_consumption,
    n_objectives=2,
    requires=[process_simulator, frac_opt],
    minimize=False,
)

Note that I added both switch_time and cycle_time as optimization variables. By adding a variable dependency, cycle_time is always kept at 4*swich_time. At some point we could consider allowing event dependencies on the cycle time but for now this is a “good enough” workaround.
Then, w_r and w_e need a special pre_processing feature. This ensures that at no point the sum of the split ratios split is larger than one. Again, this could be improved in the future by separating the assembly of actual values before setting them on FlowSheet level, but this way we always avoid inconsistent states.
For process evaluation, the fractionation optimizer is set up. It will automatically check all outlets that were specified as product_outlet in the FlowSheet and ensures that only regions with a cumulative purity higher than the one specified are accounted for.
Finally, the EluentConsumption performance indicator automatically detects Inlets that were specified as eluent_source and determines the volume that has entered the system. For it to be evaluated, it requires the process_simulator, as well as the fractionation optimizer.
Note that we suggest using multi-objective optimization.

Hope that helps. Let us know if you have further questions.

Edit: I just realized that some of these features might actually be still in beta so you would have to install the dev version of CADET-Process.

pip install git+https://github.com/fau-advanced-separations/CADET-Process.git@dev

I hope to release a new version soon™ .

By the way, we’re happy to announce our very first CADET-Workshop in the United States, hosted by our friends at RPI in Troy, NY, May 22-24, 2024. It would be great to meet you there and I’d be more than happy to do offer some more hands-on work with you on your particular problem(s). For more information on how to register, see also here!

bpimentel · April 9, 2024, 7:34pm

@j.schmoelder , thanks for the response, but I fear you may have given me more questions than answers…

So I think there may be some issues in using the CarouselBuilder vs Flowsheet, at least in this formulation of the problem. I did update to the dev version, though.

1) I’m not sure I ever specified an eluent source? Unless eluent is somehow a protected variable name

#In's
feed = Inlet(component_system, name='feed')
feed.c = [10, 10]          # M
feed.flow_rate = 2e-7    #m^3/s
eluent = Inlet(component_system, name='eluent')
eluent.c = [0, 0]          # M
eluent.flow_rate = 6e-7    #m^3/s

#Outs
raffinate = Outlet(component_system, name='raffinate')
extract = Outlet(component_system, name='extract')

2) Not all attributes seem to transfer from the FlowSheet via the CarouselBuilder wrapper.

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Cell In[25], line 1
----> 1 builder.output_states

AttributeError: 'CarouselBuilder' object has no attribute 'output_states'

So I’m not sure how to access w_r and w_e. Adding them just as names gives a different error (see #2)

3) Using builder.eluent.flow_rate since I didn’t use flow_sheet

optimization_problem.add_variable(
    'builder.eluent.flow_rate',
    lb=1e-7, ub=10e-7,
    transform='auto'
  )

raises an error

File ~\Anaconda3\envs\cadet\lib\site-packages\CADETProcess\optimization\optimizationProblem.py:347, in OptimizationProblem.add_variable(self, name, evaluation_objects, parameter_path, lb, ub, transform, indices)
    342 if parameter_path is not None and len(evaluation_objects) == 0:
    343     raise ValueError(
    344         "Cannot set parameter_path for variable without evaluation object "
    345     )
--> 347 var = OptimizationVariable(
    348     name, evaluation_objects, parameter_path,
    349     lb=lb, ub=ub, transform=transform,
    350     indices=indices,
    351 )
    353 self._variables.append(var)
    355 with warnings.catch_warnings():


CADETProcessError: Not a valid Optimization variable

Although builder.set_output_state is not an issue in the process definition section, I just can’t seem to access the value as a variable

Adding w_r and ```w_e`` give a similar issue

4) The use of the pre_processing argument in add_variable gives an unexpected keyword error

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[11], line 1
----> 1 optimization_problem.add_variable(
      2     name='w_r',
      3     lb=0, ub=1,
      4     pre_processing=lambda x: [x[0], 1 - x[0]],
      5 )
      6 optimization_problem.add_variable(
      7     name='w_e',
      8     lb=0, ub=1,
      9 )

TypeError: add_variable() got an unexpected keyword argument 'pre_processing

Adding the whole script for reference, in case I’m omitting anything important

# System Definition

from CADETProcess.processModel import ComponentSystem
from CADETProcess.processModel import Linear, Langmuir
from CADETProcess.processModel import LumpedRateModelWithoutPores
from CADETProcess.processModel import Inlet, Outlet
from CADETProcess.modelBuilder import SerialZone, ParallelZone, CarouselBuilder

from CADETProcess.simulator import Cadet
from CADETProcess.stationarity import StationarityEvaluator
from CADETProcess.stationarity import RelativeArea


from CADETProcess.processModel import Linear

from CADETProcess.optimization import OptimizationProblem
from CADETProcess.fractionation import FractionationOptimizer
from CADETProcess.performance import EluentConsumption


#Component definitions
component_system = ComponentSystem(2)

#Binding model
binding_model = Linear(component_system)
binding_model.adsorption_rate = [6, 8]
binding_model.desorption_rate = [1, 1]

#Column configuration
column = LumpedRateModelWithoutPores(component_system, name='master_column')
column.length = 0.6    # meters
column.diameter = 0.024     # meters
column.axial_dispersion = 4.4e-7     #m/s
column.total_porosity = 0.7      # v/v

column.binding_model = binding_model


# Inlet/Outlet configurations

#In's
feed = Inlet(component_system, name='feed')
feed.c = [10, 10]          # M
feed.flow_rate = 2e-7    #m^3/s
eluent = Inlet(component_system, name='eluent')
eluent.c = [0, 0]          # M
eluent.flow_rate = 6e-7    #m^3/s

#Outs
raffinate = Outlet(component_system, name='raffinate')
extract = Outlet(component_system, name='extract')

#Zones
zone_1 = SerialZone(component_system, 'zone_1', 1)
zone_2 = SerialZone(component_system, 'zone_2', 1)
zone_3 = SerialZone(component_system, 'zone_3', 1)
zone_4 = SerialZone(component_system, 'zone_4', 1)

#Carousel Build
builder = CarouselBuilder(component_system, 'smb')
builder.column = column
builder.add_unit(feed)
builder.add_unit(eluent)

builder.add_unit(raffinate)
builder.add_unit(extract)

builder.add_unit(zone_1)
builder.add_unit(zone_2)
builder.add_unit(zone_3)
builder.add_unit(zone_4)

# System Topology
builder.add_connection(eluent, zone_1)

builder.add_connection(zone_1, extract)
builder.add_connection(zone_1, zone_2)
w_e = .15
builder.set_output_state(zone_1, [w_e, 1-w_e])

builder.add_connection(zone_2, zone_3)

builder.add_connection(feed, zone_3)

builder.add_connection(zone_3, raffinate)
builder.add_connection(zone_3, zone_4)
w_r = 0.15
builder.set_output_state(zone_3, [w_r, 1-w_r])

builder.add_connection(zone_4, zone_1)

#Switch time and build
builder.switch_time = 300
process= builder.build_process()

# CSS Settings
process_simulator = Cadet()

evaluator = StationarityEvaluator()

criterion = RelativeArea()
criterion.threshold = 5e-3
evaluator.add_criterion(criterion)

process_simulator.stationarity_evaluator = evaluator
process_simulator.evaluate_stationarity = True

process_simulator.n_cycles_max = 5 #Just to speed things up
process_simulator.n_cycles_min = 1

#Optimization Definition

optimization_problem = OptimizationProblem('eluent_usage')
optimization_problem.add_evaluation_object(process)

optimization_problem.add_variable(
    'switch_time.time',
    lb=150, ub=600,
    transform='auto'
)

optimization_problem.add_variable(
    'builder.eluent.flow_rate',
    lb=1e-7, ub=10e-7,
    transform='auto'
  )

optimization_problem.add_variable(
    name='w_r',
    lb=0, ub=1,
    pre_processing=lambda x: [x[0], 1 - x[0]],
)
optimization_problem.add_variable(
    name='w_e',
    lb=0, ub=1,
)

# Setup Simulator

optimization_problem.add_evaluator(process_simulator)
# Setup Fractionator
frac_opt = FractionationOptimizer()

optimization_problem.add_evaluator(
    frac_opt,
    kwargs={'purity_required': [0.95, 0.95]}
)

# Setup Objectives
eluent_consumption = EluentConsumption()
optimization_problem.add_objective(
    eluent_consumption,
    n_objectives=2,
    requires=[process_simulator, frac_opt],
    minimize=False,
)

#Simulate Process
simulation_results = process_simulator.simulate(process)

j.schmoelder · April 10, 2024, 6:36am

Consider the CarouselBuilder as a factory pattern. It helps in configuring a Process object (including its FlowSheet) which you can then use for simulation / optimization. In future, we might add some more functionality to the CarouselBuilder s.t. it’s possible to directly use it in optimization studies but for not we don’t directly support this, unless you build your own, custom Evaluator class (which is not too difficult, actually).

You need to specify that when adding Inlets (see documentation).

flow_sheet = process.flow_sheet
flow_sheet.add_feed_inlet('feed_unit')
flow_sheet.add_eluent_inlet('eluent_unit')
flow_sheet.add_product_outlet('outlet')

alternatively, you can also specify this when adding the units to the FlowSheet in the first place.

flow_sheet.add_unit(feed_unit, feed_inlet=True)
flow_sheet.add_unit(eluent_unit, eluent_inlet=True)
flow_sheet.add_unit(raffinate_unit, product_outlet=True)
flow_sheet.add_unit(extract_unit, product_outlet=True)

Regarding your other questions, you need to make sure that the optimization variables are properly defined.

Here, it’s not the builder that controls the eluent unit but the flow_sheet attribute of the process. Use the following:

optimization_problem.add_variable(
    'flow_sheet.eluent.flow_rate',
    lb=1e-7, ub=10e-7,
    transform='auto'
  )

bpimentel:

optimization_problem.add_variable(
    name='w_r',
    lb=0, ub=1,
    pre_processing=lambda x: [x[0], 1 - x[0]],
)
optimization_problem.add_variable(
    name='w_e',
    lb=0, ub=1,
)

Note, the CarouselBuilder adds a little inlet/outlet valve unit for every zone. In practice, it’s a (very) small Cstr unit operation that shouldn’t impact the final simulation results much. Each of those automatically get added to the flow sheet with the names f"{zone.name}_inlet", and f"{zone.name}_outlet", respectively.
To optimize the output_state of the flow sheet unit, you need to specify this accordingly.

Here is the complete config that should work (at least it doesn’t throw errors for me):

# System Definition

from CADETProcess.processModel import ComponentSystem
from CADETProcess.processModel import Linear, Langmuir
from CADETProcess.processModel import LumpedRateModelWithoutPores
from CADETProcess.processModel import Inlet, Outlet
from CADETProcess.modelBuilder import SerialZone, ParallelZone, CarouselBuilder

from CADETProcess.simulator import Cadet
from CADETProcess.stationarity import StationarityEvaluator
from CADETProcess.stationarity import RelativeArea


from CADETProcess.processModel import Linear

from CADETProcess.optimization import OptimizationProblem
from CADETProcess.fractionation import FractionationOptimizer
from CADETProcess.performance import EluentConsumption


#Component definitions
component_system = ComponentSystem(2)

#Binding model
binding_model = Linear(component_system)
binding_model.adsorption_rate = [6, 8]
binding_model.desorption_rate = [1, 1]

#Column configuration
column = LumpedRateModelWithoutPores(component_system, name='master_column')
column.length = 0.6    # meters
column.diameter = 0.024     # meters
column.axial_dispersion = 4.4e-7     #m/s
column.total_porosity = 0.7      # v/v

column.binding_model = binding_model


# Inlet/Outlet configurations

#In's
feed = Inlet(component_system, name='feed')
feed.c = [10, 10]          # M
feed.flow_rate = 2e-7    #m^3/s
eluent = Inlet(component_system, name='eluent')
eluent.c = [0, 0]          # M
eluent.flow_rate = 6e-7    #m^3/s

#Outs
raffinate = Outlet(component_system, name='raffinate')
extract = Outlet(component_system, name='extract')

#Zones
zone_1 = SerialZone(component_system, 'zone_1', 1)
zone_2 = SerialZone(component_system, 'zone_2', 1)
zone_3 = SerialZone(component_system, 'zone_3', 1)
zone_4 = SerialZone(component_system, 'zone_4', 1)

#Carousel Build
builder = CarouselBuilder(component_system, 'smb')
builder.column = column
builder.add_unit(feed, feed_inlet=True)
builder.add_unit(eluent, eluent_inlet=True)

builder.add_unit(raffinate, product_outlet=True)
builder.add_unit(extract, product_outlet=True)

builder.add_unit(zone_1)
builder.add_unit(zone_2)
builder.add_unit(zone_3)
builder.add_unit(zone_4)

# System Topology
builder.add_connection(eluent, zone_1)

builder.add_connection(zone_1, extract)
builder.add_connection(zone_1, zone_2)
w_e = .15
builder.set_output_state(zone_1, [w_e, 1-w_e])

builder.add_connection(zone_2, zone_3)

builder.add_connection(feed, zone_3)

builder.add_connection(zone_3, raffinate)
builder.add_connection(zone_3, zone_4)
w_r = 0.15
builder.set_output_state(zone_3, [w_r, 1-w_r])

builder.add_connection(zone_4, zone_1)

#Switch time and build
builder.switch_time = 300
process = builder.build_process()

# CSS Settings
process_simulator = Cadet()

evaluator = StationarityEvaluator()

criterion = RelativeArea()
criterion.threshold = 5e-3
evaluator.add_criterion(criterion)

process_simulator.stationarity_evaluator = evaluator
process_simulator.evaluate_stationarity = True

process_simulator.n_cycles_max = 5 #Just to speed things up
process_simulator.n_cycles_min = 1

#Optimization Definition

optimization_problem = OptimizationProblem('eluent_usage')
optimization_problem.add_evaluation_object(process)

optimization_problem.add_variable(
    'switch_time.time',
    lb=150, ub=600,
    transform='auto'
)

optimization_problem.add_variable(
    'flow_sheet.eluent.flow_rate',
    lb=1e-7, ub=10e-7,
    transform='auto'
  )

optimization_problem.add_variable(
    name='w_r',
    parameter_path='flow_sheet.output_states.zone_3_outlet',
    lb=0, ub=1,
    pre_processing=lambda x: [x[0], 1 - x[0]],
)
optimization_problem.add_variable(
    name='w_e',
    parameter_path='flow_sheet.output_states.zone_1_outlet',
    lb=0, ub=1,
    pre_processing=lambda x: [x[0], 1 - x[0]],
)

# Setup Simulator

optimization_problem.add_evaluator(process_simulator)
# Setup Fractionator
frac_opt = FractionationOptimizer()

optimization_problem.add_evaluator(
    frac_opt,
    kwargs={'purity_required': [0.95, 0.95]}
)

# Setup Objectives
eluent_consumption = EluentConsumption()
optimization_problem.add_objective(
    eluent_consumption,
    n_objectives=2,
    requires=[process_simulator, frac_opt],
    minimize=False,
)

As a side node, we recently realized that using a (very) small Cstr model for the inlet/outlet valves massively increases the simulation runtime. We’re still trying to figure out the root cause, but we assume, it’s due to the (in theory) variable volume of the Cstr model which adds a lot of stiffness to the problem when the volume stays effectively constant. On another development branch, I implemented a hotfix that instead uses a small pfr (with just 1 axial cell) that surprisingly runs much faster. You can check it out here. Obviously, we’ll try to fix this more thoroughly in CADET-Core.

bpimentel · April 11, 2024, 5:57pm

@j.schmoelder , almost there I think - biggest issue seems to be that thr process does not contain a chromatogram to pass on to the fractionator

Setting feed/eluent/product streams at the flow_sheet level

So running your configuration as-is does throw an error on my end:

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[1], line 62
     60 builder = CarouselBuilder(component_system, 'smb')
     61 builder.column = column
---> 62 builder.add_unit(feed, feed_inlet=True)
     63 builder.add_unit(eluent, eluent_inlet=True)
     65 builder.add_unit(raffinate, product_outlet=True)

TypeError: add_unit() got an unexpected keyword argument 'feed_inlet'

feed_inlet and product_outlet configurations only seem to be settable at the Flowsheet level, which I got around by accessing builder.flowsheet

builder.flow_sheet.add_feed_inlet("feed")
builder.flow_sheet.add_eluent_inlet("eluent")
builder.flow_sheet.add_product_outlet("raffinate")
builder.flow_sheet.add_product_outlet("extract")

This seems to appropriately record the inlet and outlet assignments when I check with builder.flow_sheet.product_outlets

I then used process = builder.build_process() to put the whole thing together. Interestingly,

flow_sheet = builder.build_flowsheet()
process = Process(flow_sheet, "SMBC")

Throws a range or len type error when trying to simulate it, even though it passes all the process.check_x calls.

No chromatogram for the Fractionator implementation

At this point, everything runs, but it fails at the Fractionator implementation

fractionator = Fractionator(simulation_results)

---------
*bunch of error code*

CADETProcessError: Simulation results do not contain chromatogram

Which we can verify with

simulation_results.chromatograms
[]

Unexpected arguments*
Additionally, pre_processing is still not a valid argument for the optimization_problem.add_variable() function, not sure if thats a dev/release version issue

print(inspect.signature(optimization_problem.add_variable))
(name, evaluation_objects=-1, parameter_path=None, lb=-inf, ub=inf, transform=None, indices=None)

Optimization variable configuration
Interestingly, your configuration of the optimization variables as

optimization_problem.add_evaluation_object(process)

optimization_problem.add_variable(
    "flow_sheet.eluent.flow_rate", lb=1e-7, ub=10e-7, transform="auto"
)

seems valid, despite never having defined a flow_sheet variable. I suspect it’s actually calling the builder.flow_sheet method via the evaluation_object, which caught me off guard. But trying to change it to builder.flow_sheet.eluent.flow_rate gives a non-valid optimization variable, so it must be correct. Still haven’t actually ran the optimization, obviously

Thanks again

j.schmoelder · April 12, 2024, 5:42am

Looks like you’re not running the latest dev version of CADET-Process. You might need to uninstall and reinstall.

pip uninstall cadet-process
pip install git+https://github.com/fau-advanced-separations/CADET-Process.git@dev

Btw, v0.9.0 is also on its way. We finished the release notes yesterday and just need to review one final time before releasing. Hopefully, we can finish this today.

The flow_sheet is an attribute of the Process class. Since you used the CarouselBuilder to setup the process, it gets configured automatically.