TrustConstr optimizer (scipy) in CADET-Process

[s.zurbriggen]

I am currently using the optimization module in CADET-Process to estimate different parameters by fitting chromatography simulations onto experimental data (minimizing objective function). I was able to successfully achieve this using the U_NSGA3 optimizer from pymoo. I would like to use the TrustConstr optimizer from scipy library for the same process.

Unfortunately this results in an error message, stating that the function _minimize_trustregion_constr() got an unexpected keyword argument ‘progress_frequency’.

Please note that I also used the function check_optimization_problem(optimization_problem) to check if the problem was configured correctly and supported by the optimizer, which resulted in TRUE.

Any help would be appreciated, how I could use the TrustConstr optimizer from scipy library to solve the optimization problem.

[j.schmoelder]

Hi Samuel and welcome to the forum!

Could you please provide a script that reproduces this error?

Best

Johannes

[s.zurbriggen]

Hi Johannes,
Please see my code below.

# import tha data from csv-file
os.chdir('C:/Users/Directory')

data = pd.read_csv('Data.csv',sep=',', header=[0,])
df = pd.DataFrame(data)

voidVol = 1.148 # mL
volumetricFlowrate = 1.308  # mL/min

# CADET simulation
from CADETProcess.processModel import ComponentSystem
from CADETProcess.processModel import Inlet, Cstr, Outlet, FlowSheet, TubularReactor
from CADETProcess.processModel import Process
from CADETProcess.simulator import Cadet

# Component System
component_system = ComponentSystem(['Salt'])
Q = volumetricFlowrate *1e-6/60 # m3/s

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

cstr = Cstr(component_system, 'cstr')
cstr.c = [0]
cstr.V = 1*1e-6 # will be estimated down below
cstr.flow_rate = Q

pfr = TubularReactor(component_system, 'pfr')
pfr.length = (voidVol-0.25)*1e-6 # 0.25 = half of injection volume 0.5 mL
pfr.cross_section_area = 1
pfr.axial_dispersion = 0

outlet = Outlet(component_system, 'outlet')

# Flow Sheet
flow_sheet = FlowSheet(component_system)

flow_sheet.add_unit(inlet)
flow_sheet.add_unit(pfr)
flow_sheet.add_unit(cstr)
flow_sheet.add_unit(outlet)

flow_sheet.add_connection(inlet, pfr)
flow_sheet.add_connection(pfr, cstr)
flow_sheet.add_connection(cstr, outlet)

# Process
process = Process(flow_sheet, 'process')
process.cycle_time = 10*60

# Events & Durations
eluVol = 7.052 #mL
injVol = 5.903151 #mL
time_inj = injVol*1e-6/Q
vol_load = 0.5*1e-6 # 0.5mL

_ = process.add_event('equilibration', 'flow_sheet.inlet.c', [0], 0)
_ = process.add_event('load', 'flow_sheet.inlet.c', [1000], time_inj)
_ = process.add_event('elute', 'flow_sheet.inlet.c', [0], time_inj+(vol_load/Q))

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

time = simulation_results.solution.cstr.outlet.time #s
volume = time*Q*1e6 #mL
salt = simulation_results.solution.cstr.outlet.solution[:,0]

# Setup Optimization algorithm to estimate volume of cstr
from CADETProcess.reference import ReferenceIO
from CADETProcess.comparison import Comparator
from CADETProcess.optimization import OptimizationProblem
from CADETProcess.optimization import TrustConstr

# Setup Reference & Comparator
exp_time = df.iloc[:,0].dropna()*1e-6/Q
exp_conc = df.iloc[:,1].dropna()
reference = ReferenceIO('Salt', exp_time, exp_conc)
comparator = Comparator()
comparator.add_reference(reference)
comparator.add_difference_metric('SSE', reference, 'cstr.outlet', start=((injVol-2)*1e-6/Q), end=(11*1e-6)/Q) # start => t = V/Q
metrics = comparator.evaluate(simulation_results)
print(f"{metrics = }")

# Setup Optimization 
optimization_problem = OptimizationProblem('TN')
optimization_problem.add_evaluation_object(process)
optimization_problem.add_variable('flow_sheet.cstr.V', lb=0.01*1e-8, ub=0.5*1e-6)

simulator = Cadet()
optimization_problem.add_evaluator(simulator, cache=True)
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], frequency=1, keep_progress=True)

optimizer = TrustConstr()

print(f"{OptimizationProblem.check(optimization_problem) = }")
print(f"{optimizer.check_optimization_problem(optimization_problem) = }")

if __name__ == "__main__":
    from CADETProcess import settings
    settings.working_directory = './OPT'
    settings.temp_dir = '/dev/shm'
    print('Starting Optimization')
    
    optimization_results = optimizer.optimize(optimization_problem, x0=None)
    

Resulting in the following error:

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
C:/Users/samue/Desktop/CADET-Forum/tmp\ipykernel_2528\2395121807.py in <module>
      5     print('Starting Optimization')
      6 
----> 7     optimization_results = optimizer.optimize(optimization_problem, x0=None)
      8 

~\kcadet\CADET-Process\CADETProcess\optimization\optimizer.py in optimize(self, optimization_problem, x0, save_results, results_directory, use_checkpoint, overwrite_results_directory, exist_ok, log_level, reinit_cache, delete_cache, *args, **kwargs)
    179         start = time.time()
    180 
--> 181         self.run(self.optimization_problem, x0, *args, **kwargs)
    182 
    183         self.results.time_elapsed = time.time() - start

~\kcadet\CADET-Process\CADETProcess\optimization\scipyAdapter.py in run(self, optimization_problem, x0)
     82             warnings.filterwarnings('ignore', category=OptimizeWarning)
     83             warnings.filterwarnings('ignore', category=RuntimeWarning)
---> 84             scipy_results = optimize.minimize(
     85                 objective_function,
     86                 x0=x0,

~\anaconda3\envs\cadet\lib\site-packages\scipy\optimize\_minimize.py in minimize(fun, x0, args, method, jac, hess, hessp, bounds, constraints, tol, callback, options)
    706                               constraints, callback=callback, **options)
    707     elif meth == 'trust-constr':
--> 708         res = _minimize_trustregion_constr(fun, x0, args, jac, hess, hessp,
    709                                            bounds, constraints,
    710                                            callback=callback, **options)

TypeError: _minimize_trustregion_constr() got an unexpected keyword argument 'progress_frequency'

Thanks for you help!

Samuel

[j.schmoelder]

Thanks for the script. It does indeed seem to be a bug in CADET-Process where TrustConstr doesn’t simply ignore unknown (or rather internal) options. I have already fixed it locally and hope to make a new Hotfix release in the coming days.

[j.schmoelder]

I just released v0.7.2 which includes the fix.
Please let me know if you have further issues!

[s.zurbriggen]

Hi Johannes

Thanks for the fix! I was able to successfully use the optimizer TrustConstr() in v0.7.3 to estimate one variable (CSTR volume to simulate dispersion).

In a further step, I would like to estimate multiple variables (bed_porosity, axial_dispersion, particle_porosity) in one optimization problem using the TrustConstr() optimizer.

Unfortunately, once I want to optimize multiple variables in one optimization, I receive an error message. See the script below.

# Setup Optimization algorithm to estimate: [bed_poros, axial_disp, particle_poros]
from CADETProcess.reference import ReferenceIO
from CADETProcess.comparison import Comparator
from CADETProcess.optimization import OptimizationProblem
from CADETProcess.optimization import TrustConstr

# Setup Reference & Comparator
exp_time = df2.iloc[:,2].dropna()*1e-6/Q #mL Q=V/t <=> t=V/Q
exp_conc = df2.iloc[:,3].dropna()
reference = ReferenceIO('Salt', exp_time, exp_conc)
comparator = Comparator()
comparator.add_reference(reference)
comparator.add_difference_metric('SSE', reference, 'column.outlet', start=((injVol-2)*1e-6/Q), end=(11*1e-6)/Q) # start => t = V/Q
metrics = comparator.evaluate(simulation_results)
print(f"{metrics = }")

# Setup Optimization 
optimization_problem = OptimizationProblem('TC')
optimization_problem.add_evaluation_object(process)
optimization_problem.add_variable('flow_sheet.column.bed_porosity', lb=0.2, ub=1)
optimization_problem.add_variable('flow_sheet.column.axial_dispersion', lb=1e-8, ub=1e-5)
optimization_problem.add_variable('flow_sheet.column.particle_porosity', lb=0.2, ub=1)

simulator = Cadet()
optimization_problem.add_evaluator(simulator, cache=True)
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],
                                 frequency=1, keep_progress=True)

optimizer = TrustConstr()

if __name__ == "__main__":
    from CADETProcess import settings
    settings.working_directory = './OPT'
    settings.temp_dir = '/dev/shm'
    print('Starting Optimization')
    
    optimization_results = optimizer.optimize(optimization_problem, x0=None, save_results=True, 
                                              results_directory=None, use_checkpoint=False, 
                                              overwrite_results_directory=True, exist_ok=True, 
                                              log_level='INFO', reinit_cache=True, delete_cache=True)
    

The following error message appears:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
C:/dev/shm\ipykernel_12700\1015459542.py in <module>
      5     print('Starting Optimization')
      6 
----> 7     optimization_results = optimizer.optimize(optimization_problem, x0=None, save_results=True, 
      8                                               results_directory=None, use_checkpoint=False,
      9                                               overwrite_results_directory=True, exist_ok=True,

~\anaconda3\envs\cadet\lib\site-packages\CADETProcess\optimization\optimizer.py in optimize(self, optimization_problem, x0, save_results, results_directory, use_checkpoint, overwrite_results_directory, exist_ok, log_level, reinit_cache, delete_cache, *args, **kwargs)
    208         start = time.time()
    209 
--> 210         self.run(self.optimization_problem, x0, *args, **kwargs)
    211 
    212         self.results.time_elapsed = time.time() - start

~\anaconda3\envs\cadet\lib\site-packages\CADETProcess\optimization\scipyAdapter.py in run(self, optimization_problem, x0)
     82             warnings.filterwarnings('ignore', category=OptimizeWarning)
     83             warnings.filterwarnings('ignore', category=RuntimeWarning)
---> 84             scipy_results = optimize.minimize(
     85                 objective_function,
     86                 x0=x0,

~\anaconda3\envs\cadet\lib\site-packages\scipy\optimize\_minimize.py in minimize(fun, x0, args, method, jac, hess, hessp, bounds, constraints, tol, callback, options)
    706                               constraints, callback=callback, **options)
    707     elif meth == 'trust-constr':
--> 708         res = _minimize_trustregion_constr(fun, x0, args, jac, hess, hessp,
    709                                            bounds, constraints,
    710                                            callback=callback, **options)

~\anaconda3\envs\cadet\lib\site-packages\scipy\optimize\_trustregion_constr\minimize_trustregion_constr.py in _minimize_trustregion_constr(fun, x0, args, grad, hess, hessp, bounds, constraints, xtol, gtol, barrier_tol, sparse_jacobian, callback, maxiter, verbose, finite_diff_rel_step, initial_constr_penalty, initial_tr_radius, initial_barrier_parameter, initial_barrier_tolerance, factorization_method, disp)
    507 
    508     elif method == 'tr_interior_point':
--> 509         _, result = tr_interior_point(
    510             objective.fun, objective.grad, lagrangian_hess,
    511             n_vars, canonical.n_ineq, canonical.n_eq,

~\anaconda3\envs\cadet\lib\site-packages\scipy\optimize\_trustregion_constr\tr_interior_point.py in tr_interior_point(fun, grad, lagr_hess, n_vars, n_ineq, n_eq, constr, jac, x0, fun0, grad0, constr_ineq0, jac_ineq0, constr_eq0, jac_eq0, stop_criteria, enforce_feasibility, xtol, state, initial_barrier_parameter, initial_tolerance, initial_penalty, initial_trust_radius, factorization_method)
    319     while True:
    320         # Solve SQP subproblem
--> 321         z, state = equality_constrained_sqp(
    322             subprob.function_and_constraints,
    323             subprob.gradient_and_jacobian,

~\anaconda3\envs\cadet\lib\site-packages\scipy\optimize\_trustregion_constr\equality_constrained_sqp.py in equality_constrained_sqp(fun_and_constr, grad_and_jac, lagr_hess, x0, fun0, grad0, constr0, jac0, stop_criteria, state, initial_penalty, initial_trust_radius, factorization_method, trust_lb, trust_ub, scaling)
     91 
     92     last_iteration_failed = False
---> 93     while not stop_criteria(state, x, last_iteration_failed,
     94                             optimality, constr_violation,
     95                             trust_radius, penalty, cg_info):

~\anaconda3\envs\cadet\lib\site-packages\scipy\optimize\_trustregion_constr\tr_interior_point.py in stop_criteria(self, state, z, last_iteration_failed, optimality, constr_violation, trust_radius, penalty, cg_info)
    249         """
    250         x = self.get_variables(z)
--> 251         if self.global_stop_criteria(state, x,
    252                                      last_iteration_failed,
    253                                      trust_radius, penalty,

~\anaconda3\envs\cadet\lib\site-packages\scipy\optimize\_trustregion_constr\minimize_trustregion_constr.py in stop_criteria(state, x, last_iteration_failed, tr_radius, constr_penalty, cg_info, barrier_parameter, barrier_tolerance)
    467             state.status = None
    468             state.niter = state.nit  # Alias for callback (backward compatibility)
--> 469             if callback is not None and callback(np.copy(state.x), state):
    470                 state.status = 3
    471             elif state.optimality < gtol and state.constr_violation < gtol:

~\anaconda3\envs\cadet\lib\site-packages\CADETProcess\optimization\scipyAdapter.py in callback_function(x, state)
     72             cv = optimization_problem.evaluate_nonlinear_constraints_violation(x)
     73 
---> 74             self.run_post_evaluation_processing(x, f, g, cv, self.n_evals)
     75 
     76             return False

~\anaconda3\envs\cadet\lib\site-packages\CADETProcess\optimization\optimizer.py in run_post_evaluation_processing(self, x, f, g, cv, current_evaluation, x_opt)
    380             self.results.update_pareto(pareto_front)
    381 
--> 382         self._run_post_processing(current_evaluation)
    383 
    384         self.logger.info(

~\anaconda3\envs\cadet\lib\site-packages\CADETProcess\optimization\optimizer.py in _run_post_processing(self, current_iteration)
    308 
    309         if current_iteration % self.progress_frequency == 0:
--> 310             self.results.plot_figures(show=False)
    311 
    312         for callback in self.optimization_problem.callbacks:

~\anaconda3\envs\cadet\lib\site-packages\CADETProcess\optimization\results.py in plot_figures(self, show)
    381             )
    382             if self.optimization_problem.n_variables > 1:
--> 383                 self.plot_corner(
    384                     show=show, plot_directory=self.plot_directory
    385                 )

~\anaconda3\envs\cadet\lib\site-packages\CADETProcess\optimization\results.py in plot_corner(self, *args, **kwargs)
    536         """
    537         try:
--> 538             self.population_all.plot_corner(*args, **kwargs)
    539         except AssertionError:
    540             pass

~\anaconda3\envs\cadet\lib\site-packages\CADETProcess\optimization\population.py in plot_corner(self, untransformed, show, plot_directory)
    573         labels = [label for label in labels if label not in singular_labels]
    574 
--> 575         fig = corner.corner(
    576             x,
    577             labels=labels,

~\anaconda3\envs\cadet\lib\site-packages\corner\corner.py in corner(data, bins, range, weights, color, hist_bin_factor, smooth, smooth1d, labels, label_kwargs, titles, show_titles, title_fmt, title_kwargs, truths, truth_color, scale_hist, quantiles, verbose, fig, max_n_ticks, top_ticks, use_math_text, reverse, labelpad, hist_kwargs, group, var_names, filter_vars, coords, divergences, divergences_kwargs, labeller, **hist2d_kwargs)
    254         )
    255 
--> 256     return arviz_corner(
    257         data,
    258         bins=bins,

~\anaconda3\envs\cadet\lib\site-packages\corner\arviz_corner.py in arviz_corner(data, bins, range, weights, color, hist_bin_factor, smooth, smooth1d, labels, label_kwargs, titles, show_titles, title_fmt, title_kwargs, truths, truth_color, scale_hist, quantiles, verbose, fig, max_n_ticks, top_ticks, use_math_text, reverse, labelpad, hist_kwargs, group, var_names, filter_vars, coords, divergences, divergences_kwargs, labeller, **hist2d_kwargs)
    133 
    134     # Coerce the samples into the expected format
--> 135     samples = np.stack([x[-1].flatten() for x in plotters], axis=-1)
    136     fig = corner_impl(
    137         samples,

~\anaconda3\envs\cadet\lib\site-packages\numpy\core\overrides.py in stack(*args, **kwargs)

~\anaconda3\envs\cadet\lib\site-packages\numpy\core\shape_base.py in stack(arrays, axis, out)
    420     arrays = [asanyarray(arr) for arr in arrays]
    421     if not arrays:
--> 422         raise ValueError('need at least one array to stack')
    423 
    424     shapes = {arr.shape for arr in arrays}

ValueError: need at least one array to stack

How can I optimze the three variables in one optimization using the TrustConstr() optimizer?

Your help is greatly appreciated!

[alters]

In practice, those three variables should be measured experimentally (e.g., tracer experiments for porosity) or using correlations for axial dispersion (e.g., Rastegar and Gu 2017 “Empirical correlations for axial dispersion coefficient and Pecletnumber in fixed-bed columns”).

Intraparticle porosity (particle porosity) values have also been reported in the literature for various stationary phases and interstitial porosity (bed porosity). I have never seen interstitial porosity outside of the range of 0.3 - 0.45. In fact, interstitial porosity is bounded by the theoretical limits of random sphere packing (0.26-0.48). It is a reasonable assumption to use 0.4 for interstitial porosity since this is a pretty typical value for packed columns. Intraparticle porosity has a much bigger range and is dependent on the tracer used, but for fully pore-penetrating tracers (e.g., salt, acetone, or glucose) it typically ranges between 0.6 and 0.9.

I don’t recommend fitting because the value may not be determinable and may give you something nonphysical, depending on what bounds you use.

[j.schmoelder]

Hey Samuel,

unfortunately, it’s not possible for me to reproduce your issue. Please provide all required data and scripts in a minimal reproducible example.

Best

Johannes