Linear equality constraint - Cadet Process

CADET Troubleshooting
1

Hello,

I am having problems implementing a linear equality constraint. I am trying to implement that the axial dispersion of two tubings, which are both optimisation variables, should be identical.
I have tried a few things for this. The one I thought would work most likely is this:

optimization_problem.add_variable(
    name='axial_dispersion_column', parameter_path='flow_sheet.column.axial_dispersion',
    lb=1e-10, ub=0.1,
    transform='auto'
)

optimization_problem.add_variable(
    name='axial_dispersion_tubing', parameter_path='flow_sheet.tubing1.axial_dispersion',
    lb=1e-10, ub=0.1,
    transform='auto'
)

optimization_problem.add_variable(
    name='axial_dispersion_tubing2', parameter_path='flow_sheet.tubing2.axial_dispersion',
    lb=1e-10, ub=0.1,
    transform='auto'
)


optimization_problem.add_linear_equality_constraint(['axial_dispersion_tubing', 'axial_dispersion_tubing2'], [1, -1], 0)

However, it does not work, which can also be seen from the fact that for the following method:

optimization_problem.check_linear_equality_constraints([4e-08, 8e-06, 8e-06])

I get the following output: False

Where is the error in my code? Many thanks in advance.

Stefanos

2

Brief update:
I have found the cause. It works with the following line of code:
optimization_problem.add_linear_equality_constraint(['axial_dispersion_tubing', 'axial_dispersion_tubing2'], lhs=[1, -1], beq=0, eps=1e-9)

However, I now get an error during optimization:

CADETProcessError                         Traceback (most recent call last)
Cell In[14], line 2
      1 if __name__ == '__main__':
----> 2     optimization_results = optimizer.optimize(
      3         optimization_problem, x0,
      4         use_checkpoint=True
      5     )

File ~\anaconda3\envs\MasterNew\lib\site-packages\CADETProcess\optimization\optimizer.py:260, in OptimizerBase.optimize(self, optimization_problem, x0, save_results, results_directory, use_checkpoint, overwrite_results_directory, exist_ok, log_level, reinit_cache, delete_cache, *args, **kwargs)
    256 plt.switch_backend('agg')
    258 start = time.time()
--> 260 self.run(self.optimization_problem, x0, *args, **kwargs)
    262 self.results.time_elapsed = time.time() - start
    264 if delete_cache:

File ~\anaconda3\envs\MasterNew\lib\site-packages\CADETProcess\optimization\pymooAdapter.py:92, in PymooInterface.run(self, optimization_problem, x0)
     90 if len(pop) < pop_size:
     91     n_remaining = pop_size - len(pop)
---> 92     remaining = optimization_problem.create_initial_values(
     93         n_remaining, method='chebyshev', seed=self.seed,
     94         burn_in=burn_in
     95     )
     96     pop = np.vstack((pop, remaining))
     97 elif len(pop) > pop_size:

File ~\anaconda3\envs\MasterNew\lib\site-packages\CADETProcess\optimization\optimizationProblem.py:2780, in OptimizationProblem.create_initial_values(self, n_samples, method, seed, burn_in)
   2778 while len(values) < n_samples:
   2779     if counter > burn_in:
-> 2780         raise CADETProcessError(
   2781             "Cannot find invididuals that fulfill constraints."
   2782         )
   2784     counter += 1
   2785     i = rng.integers(0, burn_in)

CADETProcessError: Cannot find invididuals that fulfill constraints.

I think it’s best if I share the whole code with the following experimental data:
NaCl_Normalized.csv (63.0 KB)

from CADETProcess.processModel import ComponentSystem
from CADETProcess.processModel import StericMassAction
from CADETProcess.processModel import Inlet, GeneralRateModel, Outlet,TubularReactor, Cstr,LumpedRateModelWithPores
from CADETProcess.processModel import FlowSheet
from CADETProcess.processModel import Process
from CADETProcess.simulator import Cadet
import numpy as np

Q_ml_min = 0.5  # ml/min
Q_m3_s = Q_ml_min/(60*1e6) #m3/s

V_tracer = 5e-8  

component_system = ComponentSystem()
component_system.add_component('Salt')


# Unit Operations
inlet = Inlet(component_system, name='inlet')
inlet.flow_rate = Q_m3_s


column = LumpedRateModelWithPores(component_system, 'column')

column.length = 0.1
column.diameter = 0.0077
column.bed_porosity = 0.19 
column.particle_radius = 3.4e-5 
column.particle_porosity = 0.8333
column.axial_dispersion =4.41924993e-08
column.film_diffusion = [1.829e-4]


# Tubing1
tubing1 = TubularReactor(component_system, 'tubing1')
tubing1.length = 0.534
tubing1.cross_section_area = 4.418e-7
tubing1.axial_dispersion = 8.97920728e-04 


tubing2 = TubularReactor(component_system, 'tubing2')
tubing2.length = 0.534 
tubing2.cross_section_area = 4.418e-7 
tubing2.axial_dispersion =   8.97920728e-04 


uv_detector = Cstr(component_system, name='uv_detector')
uv_detector.V = 2e-9 
uv_detector.c = [0]


outlet = Outlet(component_system, name='outlet')

# Flow Sheet
flow_sheet = FlowSheet(component_system)


flow_sheet.add_unit(inlet)
flow_sheet.add_unit(column)
flow_sheet.add_unit(outlet)
flow_sheet.add_unit(tubing1)
flow_sheet.add_unit(tubing2)
flow_sheet.add_unit(uv_detector)

flow_sheet.add_connection(inlet, tubing1)
flow_sheet.add_connection(tubing1, column)
flow_sheet.add_connection(column, tubing2)
flow_sheet.add_connection(tubing2, uv_detector)
flow_sheet.add_connection(uv_detector, outlet)


# Process
process = Process(flow_sheet, 'lwe')
process.cycle_time = 1200

_ = process.add_event('Pre-Wash', 'flow_sheet.inlet.c', [0],0)
_ = process.add_event('Injection', 'flow_sheet.inlet.c', [500], 60)
_ = process.add_event('Wash', 'flow_sheet.inlet.c', [0],(V_tracer/Q_m3_s+60))

process_simulator = Cadet()

simulation_results = process_simulator.simulate(process)
fig, ax = simulation_results.solution.outlet.outlet.plot()

from CADETProcess.reference import ReferenceIO
import numpy as np

data_uv_2 = np.loadtxt('NaCl_Normalized.csv', delimiter=',')
time_uv_2 = data_uv_2[:, 0] * 60
uv_2 = data_uv_2[:, 1]

reference = ReferenceIO('NaCl experiment', time_uv_2, uv_2) 
_ = reference.plot()

from CADETProcess.comparison import Comparator
comparator = Comparator()

comparator.add_reference(reference)
comparator.add_difference_metric(
        'PeakPosition', reference,
        'outlet.outlet'
    )
comparator.add_difference_metric(
        'PeakHeight', reference,
        'outlet.outlet'
    )
from CADETProcess.simulator import Cadet
simulator = Cadet()
simulation_results = simulator.simulate(process)

metrics = comparator.evaluate(simulation_results)
print(metrics)
plot= comparator.plot_comparison(simulation_results) 

from CADETProcess.optimization import OptimizationProblem
optimization_problem = OptimizationProblem('AxialDispersion_PositionHeight')

optimization_problem.add_evaluation_object(process)



optimization_problem.add_variable(
    name='axial_dispersion_column', parameter_path='flow_sheet.column.axial_dispersion',
    lb=1e-10, ub=0.1,
    transform='auto'
)

optimization_problem.add_variable(
    name='axial_dispersion_tubing', parameter_path='flow_sheet.tubing1.axial_dispersion',
    lb=1e-9, ub=0.1,
    transform='auto'
)

optimization_problem.add_variable(
    name='axial_dispersion_tubing2', parameter_path='flow_sheet.tubing2.axial_dispersion',
    lb=1e-9, ub=0.1,
    transform='auto'
)

optimization_problem.add_linear_equality_constraint(['axial_dispersion_tubing', 'axial_dispersion_tubing2'], lhs=[1, -1], beq=0, eps=1e-9)
#optimization_problem.add_linear_equality_constraint(['axial_dispersion_tubing', 'axial_dispersion_tubing2'], lhs=[1e5, -1e5], beq=0, eps=1e-4) # I have also tried this


optimization_problem.add_evaluator(simulator)

optimization_problem.add_objective(
    comparator,
    n_objectives=comparator.n_metrics,
    requires=[simulator]
)

def callback(simulation_results, individual, evaluation_object, callbacks_dir='./'):
    comparator.plot_comparison(
        simulation_results,
        file_name=f'{callbacks_dir}/{individual.id}_{evaluation_object}_comparison.png',
        show=False
    )


optimization_problem.add_callback(callback, requires=[simulator])

x0=[[5e-08, 5e-06,5e-06]]
#x0=optimization_problem.create_initial_values(method='chebyshev')

from CADETProcess.optimization import U_NSGA3
optimizer = U_NSGA3()

optimizer.n_cores = 4 
optimizer.pop_size = 50
optimizer.n_max_gen = 5

if __name__ == '__main__':
    optimization_results = optimizer.optimize(
        optimization_problem, x0,
        use_checkpoint=False
    )

Epsilon must be chosen very low, as the axial dispersion is very small. As shown above, I have also tried multiplication with a pre-factor and other possibilities (e.g. with initial values generated by the appropriate method).

Can anyone help me? :slight_smile:

3

I’ve found the source, will work on a fix later today or tomorrow.

4

Hey,

so we would recommend to implement these types of parameter dependencies using “optimization_problem.add_variable_dependency”.

optimization_problem.add_variable_dependency(dependent_variable="axial_dispersion_tubing2",
                                             independent_variables="axial_dispersion_tubing",
                                             transform=lambda x: x)

This should work:

from CADETProcess.processModel import ComponentSystem
from CADETProcess.processModel import StericMassAction
from CADETProcess.processModel import Inlet, GeneralRateModel, Outlet, TubularReactor, Cstr, LumpedRateModelWithPores
from CADETProcess.processModel import FlowSheet
from CADETProcess.processModel import Process
from CADETProcess.simulator import Cadet
import numpy as np

Q_ml_min = 0.5  # ml/min
Q_m3_s = Q_ml_min / (60 * 1e6)  # m3/s

V_tracer = 5e-8

component_system = ComponentSystem()
component_system.add_component('Salt')

# Unit Operations
inlet = Inlet(component_system, name='inlet')
inlet.flow_rate = Q_m3_s

column = LumpedRateModelWithPores(component_system, 'column')

column.length = 0.1
column.diameter = 0.0077
column.bed_porosity = 0.19
column.particle_radius = 3.4e-5
column.particle_porosity = 0.8333
column.axial_dispersion = 4.41924993e-08
column.film_diffusion = [1.829e-4]

# Tubing1
tubing1 = TubularReactor(component_system, 'tubing1')
tubing1.length = 0.534
tubing1.cross_section_area = 4.418e-7
tubing1.axial_dispersion = 8.97920728e-04

tubing2 = TubularReactor(component_system, 'tubing2')
tubing2.length = 0.534
tubing2.cross_section_area = 4.418e-7
tubing2.axial_dispersion = 8.97920728e-04

uv_detector = Cstr(component_system, name='uv_detector')
uv_detector.V = 2e-9
uv_detector.c = [0]

outlet = Outlet(component_system, name='outlet')

# Flow Sheet
flow_sheet = FlowSheet(component_system)

flow_sheet.add_unit(inlet)
flow_sheet.add_unit(column)
flow_sheet.add_unit(outlet)
flow_sheet.add_unit(tubing1)
flow_sheet.add_unit(tubing2)
flow_sheet.add_unit(uv_detector)

flow_sheet.add_connection(inlet, tubing1)
flow_sheet.add_connection(tubing1, column)
flow_sheet.add_connection(column, tubing2)
flow_sheet.add_connection(tubing2, uv_detector)
flow_sheet.add_connection(uv_detector, outlet)

# Process
process = Process(flow_sheet, 'lwe')
process.cycle_time = 1200

_ = process.add_event('Pre-Wash', 'flow_sheet.inlet.c', [0], 0)
_ = process.add_event('Injection', 'flow_sheet.inlet.c', [500], 60)
_ = process.add_event('Wash', 'flow_sheet.inlet.c', [0], (V_tracer / Q_m3_s + 60))

process_simulator = Cadet()
process_simulator.time_resolution = 10

simulation_results = process_simulator.simulate(process)
fig, ax = simulation_results.solution.outlet.outlet.plot()

from CADETProcess.reference import ReferenceIO
import numpy as np

data_uv_2 = np.loadtxt('NaCl_Normalized.csv', delimiter=',')
time_uv_2 = data_uv_2[:, 0] * 60
uv_2 = data_uv_2[:, 1]

reference = ReferenceIO('NaCl experiment', time_uv_2, uv_2)
_ = reference.plot()

from CADETProcess.comparison import Comparator

comparator = Comparator()

comparator.add_reference(reference)
comparator.add_difference_metric(
    'PeakPosition', reference,
    'outlet.outlet'
)
comparator.add_difference_metric(
    'PeakHeight', reference,
    'outlet.outlet'
)
from CADETProcess.simulator import Cadet

simulator = Cadet()
simulation_results = simulator.simulate(process)

metrics = comparator.evaluate(simulation_results)
print(metrics)
plot = comparator.plot_comparison(simulation_results)

from CADETProcess.optimization import OptimizationProblem

optimization_problem = OptimizationProblem('AxialDispersion_PositionHeight')

optimization_problem.add_evaluation_object(process)

optimization_problem.add_variable(
    name='axial_dispersion_column', parameter_path='flow_sheet.column.axial_dispersion',
    lb=1e-10, ub=0.1,
    transform='auto'
)

optimization_problem.add_variable(
    name='axial_dispersion_tubing', parameter_path='flow_sheet.tubing1.axial_dispersion',
    lb=1e-9, ub=0.1,
    transform='auto'
)

optimization_problem.add_variable(
    name='axial_dispersion_tubing2', parameter_path='flow_sheet.tubing2.axial_dispersion',
    lb=1e-9, ub=0.1,
    transform='auto'
)

optimization_problem.add_variable_dependency(dependent_variable="axial_dispersion_tubing2",
                                             independent_variables="axial_dispersion_tubing",
                                             transform=lambda x: x)

optimization_problem.add_evaluator(simulator)

optimization_problem.add_objective(
    comparator,
    n_objectives=comparator.n_metrics,
    requires=[simulator]
)


def callback(simulation_results, individual, evaluation_object, callbacks_dir='./'):
    comparator.plot_comparison(
        simulation_results,
        file_name=f'{callbacks_dir}/{individual.id}_{evaluation_object}_comparison.png',
        show=False
    )


optimization_problem.add_callback(callback, requires=[simulator])

x0 = [[5e-08, 5e-06, 5e-06]]
# x0=optimization_problem.create_initial_values(method='chebyshev')

from CADETProcess.optimization import U_NSGA3

optimizer = U_NSGA3()

optimizer.n_cores = 4
optimizer.pop_size = 4
optimizer.n_max_gen = 2

import matplotlib.pyplot as plt

plt.close("all")

if __name__ == '__main__':
    optimization_results = optimizer.optimize(
        optimization_problem, x0,
        use_checkpoint=False
    )

We did discover another bug in the code that deals with variable dependencies, so you’ll have to update CADET Process again to the current dev version before using that code.

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