TrustConstr optimizer (scipy) in CADET-Process

CADET Troubleshooting
1

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.

2

Hi Samuel and welcome to the forum!

Could you please provide a script that reproduces this error?

Best

Johannes

3

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

4

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.

5

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

1 Like
6

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!

7

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.

1 Like
8

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

9

Hey Johannes

Thanks for the CADET Workshop last week. It was very helpful. Please find below the minimal reproducible example.

RefExp.csv (102.8 KB)

import numpy as np
from scipy.interpolate import interp1d
from scipy.signal import find_peaks

import matplotlib.pyplot as plt 
import pandas as pd
import os

from CADETProcess.processModel import (
    ComponentSystem, Inlet, Cstr, TubularReactor, Outlet, FlowSheet, Process)

from CADETProcess.simulator import Cadet
from CADETProcess.reference import ReferenceIO
from CADETProcess.comparison import Comparator
from CADETProcess.optimization import OptimizationProblem, TrustConstr

aceton_1 = pd.read_csv("RefExp.csv",delimiter=",", header=0)

#Variables
Q_mLmin = 0.982
Q_m3s = Q_mLmin/(60*1e6)
V_tracer = 0.5*1e-6

#Component System 
component_system = ComponentSystem(['Aceton'])

#Unit Operations
inlet_aceton = Inlet(component_system, name='inlet_aceton')
inlet_aceton.c = [10*1e3/58.08]

inlet_water = Inlet(component_system, name='inlet_water')
inlet_water.c = [0]

mixer = Cstr(component_system, name='mixer')
mixer.V = 1.37910333e-07
mixer.c = [0]

tubing = TubularReactor(component_system, name='tubing')
tubing.length = 0.898*1e-6
tubing.cross_section_area = 1
tubing.axial_dispersion = 0
tubing.c = [0]

uv_det = Cstr(component_system, name='uv_det')
uv_det.V = 0.03*1e-6
uv_det.c = [0]

cond_det = Cstr(component_system, name='cond_det')
cond_det.V = 0.03*1e-6
cond_det.c = [0]

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

#FlowSheet
flowsheet = FlowSheet(component_system)
flowsheet.add_unit(inlet_aceton)
flowsheet.add_unit(inlet_water)
flowsheet.add_unit(mixer)
flowsheet.add_unit(tubing)
flowsheet.add_unit(uv_det)
flowsheet.add_unit(cond_det)
flowsheet.add_unit(outlet)

flowsheet.add_connection(inlet_aceton, mixer)
flowsheet.add_connection(inlet_water, mixer)
flowsheet.add_connection(mixer, tubing)
flowsheet.add_connection(tubing, uv_det)
flowsheet.add_connection(uv_det, cond_det)
flowsheet.add_connection(cond_det, outlet)

#Process
process = Process(flowsheet, name='process')
process.cycle_time = 1000

start_equil_water = 10*1e-6/((5/0.982)*Q_m3s)
start_salt_on = start_equil_water + 4.91*1e-6/(Q_m3s)
start_water_flush = start_salt_on + V_tracer/Q_m3s

process.add_event('pre-equil_water', 'flow_sheet.inlet_water.flow_rate', state=(5/0.982)*Q_m3s, time=0)
process.add_event('equil_water', 'flow_sheet.inlet_water.flow_rate', state=Q_m3s, time=start_equil_water)
process.add_event('salt_on', 'flow_sheet.inlet_aceton.flow_rate', state=Q_m3s, time=start_salt_on)
process.add_event('water_flush', 'flow_sheet.inlet_water.flow_rate', state=Q_m3s, time=start_water_flush) 

process.add_event('salt_off', 'flow_sheet.inlet_aceton.flow_rate', state=0, dependencies=['water_flush'])
process.add_event('water_off', 'flow_sheet.inlet_water.flow_rate', state=0, dependencies=['salt_on'])

#Simulation
process_simulator = Cadet()

simulation_results = process_simulator.simulate(process)

>
>Reference and Comparator
reference_aceton_1 = ReferenceIO('Aceton_1', aceton_1.iloc[:,0]*60, aceton_1.iloc[:,1])

comparator = Comparator()
comparator.add_reference(reference_aceton_1)

comparator.add_difference_metric('SSE', reference_aceton_1, 'uv_det.outlet', start=7*60, end=10.5*60)
comparator.evaluate(simulation_results)

comparator.plot_comparison(simulation_results)
>
>
dir = 'G:/XX'
os.chdir(dir)
cwd = os.getcwd(); print(cwd)

optimization_problem = OptimizationProblem('PeakNoCol')
optimization_problem.add_evaluation_object(process)
optimization_problem.add_variable(name='flow_sheet.mixer.V', lb=0.01*1e-6, ub=1.5*1e-6)
optimization_problem.add_variable(name='flow_sheet.tubing.axial_dispersion', lb=1e-19, ub=1e-10)

optimization_problem.add_evaluator(process_simulator, cache=True)
optimization_problem.add_objective(comparator, n_objectives=comparator.n_metrics, requires=[process_simulator])

optimizer = TrustConstr()
optimizer.n_cores = 1
>
 optimization_results = optimizer.optimize(optimization_problem, x0=[0.5e-06, 1e-16], save_results=True, results_directory=dir, use_checkpoint=False, overwrite_results_directory=False, exist_ok=True, log_level='INFO', reinit_cache=False, delete_cache=True)

Resulting in the follwoing error:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[16], line 1
----> 1 optimization_results = optimizer.optimize(optimization_problem, x0=[0.5e-06, 1e-16], save_results=True, results_directory=dir, use_checkpoint=False,
      2                                          overwrite_results_directory=False, exist_ok=True, log_level='INFO', reinit_cache=False, delete_cache=True)

File C:\ProgramData\miniconda3\envs\cadet_env\lib\site-packages\CADETProcess\optimization\optimizer.py:210, 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)
    206 plt.switch_backend('agg')
    208 start = time.time()
--> 210 self.run(self.optimization_problem, x0, *args, **kwargs)
    212 self.results.time_elapsed = time.time() - start
    214 if delete_cache:

File C:\ProgramData\miniconda3\envs\cadet_env\lib\site-packages\CADETProcess\optimization\scipyAdapter.py:84, in SciPyInterface.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,
     87         method=str(self),
     88         tol=self.tol,
     89         jac=self.jac,
     90         constraints=self.get_constraint_objects(optimization_problem),
     91         bounds=self.get_bounds(optimization_problem),
     92         options=self.specific_options,
     93         callback=callback_function,
     94     )
     96 self.results.exit_flag = scipy_results.status
     97 self.results.exit_message = scipy_results.message

File C:\ProgramData\miniconda3\envs\cadet_env\lib\site-packages\scipy\optimize\_minimize.py:722, in minimize(fun, x0, args, method, jac, hess, hessp, bounds, constraints, tol, callback, options)
    719     res = _minimize_slsqp(fun, x0, args, jac, bounds,
    720                           constraints, callback=callback, **options)
    721 elif meth == 'trust-constr':
--> 722     res = _minimize_trustregion_constr(fun, x0, args, jac, hess, hessp,
    723                                        bounds, constraints,
    724                                        callback=callback, **options)
    725 elif meth == 'dogleg':
    726     res = _minimize_dogleg(fun, x0, args, jac, hess,
    727                            callback=callback, **options)

File C:\ProgramData\miniconda3\envs\cadet_env\lib\site-packages\scipy\optimize\_trustregion_constr\minimize_trustregion_constr.py:528, 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)
    519     _, result = equality_constrained_sqp(
    520         fun_and_constr, grad_and_jac, lagrangian_hess,
    521         x0, objective.f, objective.g,
   (...)
    524         initial_constr_penalty, initial_tr_radius,
    525         factorization_method)
    527 elif method == 'tr_interior_point':
--> 528     _, result = tr_interior_point(
    529         objective.fun, objective.grad, lagrangian_hess,
    530         n_vars, canonical.n_ineq, canonical.n_eq,
    531         canonical.fun, canonical.jac,
    532         x0, objective.f, objective.g,
    533         c_ineq0, J_ineq0, c_eq0, J_eq0,
    534         stop_criteria,
    535         canonical.keep_feasible,
    536         xtol, state, initial_barrier_parameter,
    537         initial_barrier_tolerance,
    538         initial_constr_penalty, initial_tr_radius,
    539         factorization_method)
    541 # Status 3 occurs when the callback function requests termination,
    542 # this is assumed to not be a success.
    543 result.success = True if result.status in (1, 2) else False

File C:\ProgramData\miniconda3\envs\cadet_env\lib\site-packages\scipy\optimize\_trustregion_constr\tr_interior_point.py:321, 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)
    318 # Solves a sequence of barrier problems
    319 while True:
    320     # Solve SQP subproblem
--> 321     z, state = equality_constrained_sqp(
    322         subprob.function_and_constraints,
    323         subprob.gradient_and_jacobian,
    324         subprob.lagrangian_hessian,
    325         z, fun0_subprob, grad0_subprob,
    326         constr0_subprob, jac0_subprob, subprob.stop_criteria,
    327         state, initial_penalty, trust_radius,
    328         factorization_method, trust_lb, trust_ub, subprob.scaling)
    329     if subprob.terminate:
    330         break

File C:\ProgramData\miniconda3\envs\cadet_env\lib\site-packages\scipy\optimize\_trustregion_constr\equality_constrained_sqp.py:93, 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)
     89 cg_info = {'niter': 0, 'stop_cond': 0,
     90            'hits_boundary': False}
     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):
     96     # Normal Step - `dn`
     97     # minimize 1/2*||A dn + b||^2
     98     # subject to:
     99     # ||dn|| <= TR_FACTOR * trust_radius
    100     # BOX_FACTOR * lb <= dn <= BOX_FACTOR * ub.
    101     dn = modified_dogleg(A, Y, b,
    102                          TR_FACTOR*trust_radius,
    103                          BOX_FACTOR*trust_lb,
    104                          BOX_FACTOR*trust_ub)
    106     # Tangential Step - `dt`
    107     # Solve the QP problem:
    108     # minimize 1/2 dt.T H dt + dt.T (H dn + c)
   (...)
    111     # ||dt|| <= sqrt(trust_radius**2 - ||dn||**2)
    112     # lb - dn <= dt <= ub - dn

File C:\ProgramData\miniconda3\envs\cadet_env\lib\site-packages\scipy\optimize\_trustregion_constr\tr_interior_point.py:251, in BarrierSubproblem.stop_criteria(self, state, z, last_iteration_failed, optimality, constr_violation, trust_radius, penalty, cg_info)
    246 """Stop criteria to the barrier problem.
    247 The criteria here proposed is similar to formula (2.3)
    248 from [1]_, p.879.
    249 """
    250 x = self.get_variables(z)
--> 251 if self.global_stop_criteria(state, x,
    252                              last_iteration_failed,
    253                              trust_radius, penalty,
    254                              cg_info,
    255                              self.barrier_parameter,
    256                              self.tolerance):
    257     self.terminate = True
    258     return True

File C:\ProgramData\miniconda3\envs\cadet_env\lib\site-packages\scipy\optimize\_trustregion_constr\minimize_trustregion_constr.py:484, in _minimize_trustregion_constr.<locals>.stop_criteria(state, x, last_iteration_failed, tr_radius, constr_penalty, cg_info, barrier_parameter, barrier_tolerance)
    482 callback_stop = False
    483 try:
--> 484     callback_stop = callback(state)
    485 except StopIteration:
    486     callback_stop = True

File C:\ProgramData\miniconda3\envs\cadet_env\lib\site-packages\scipy\optimize\_optimize.py:147, in _wrap_callback.<locals>.wrapped_callback(res)
    146 def wrapped_callback(res):
--> 147     return callback(np.copy(res.x), res)

File C:\ProgramData\miniconda3\envs\cadet_env\lib\site-packages\CADETProcess\optimization\scipyAdapter.py:74, in SciPyInterface.run.<locals>.callback_function(x, state)
     71 g = optimization_problem.evaluate_nonlinear_constraints(x)
     72 cv = optimization_problem.evaluate_nonlinear_constraints_violation(x)
---> 74 self.run_post_evaluation_processing(x, f, g, cv, self.n_evals)
     76 return False

File C:\ProgramData\miniconda3\envs\cadet_env\lib\site-packages\CADETProcess\optimization\optimizer.py:382, in OptimizerBase.run_post_evaluation_processing(self, x, f, g, cv, current_evaluation, x_opt)
    378     pareto_front.add_individual(ind)
    380     self.results.update_pareto(pareto_front)
--> 382 self._run_post_processing(current_evaluation)
    384 self.logger.info(
    385     f'Finished Evaluation {current_evaluation}.'
    386 )
    387 for ind in self.results.pareto_front:

File C:\ProgramData\miniconda3\envs\cadet_env\lib\site-packages\CADETProcess\optimization\optimizer.py:310, in OptimizerBase._run_post_processing(self, current_iteration)
    307     self.results.update_meta(meta_front)
    309 if current_iteration % self.progress_frequency == 0:
--> 310     self.results.plot_figures(show=False)
    312 for callback in self.optimization_problem.callbacks:
    313     if self.optimization_problem.n_callbacks > 1:

File C:\ProgramData\miniconda3\envs\cadet_env\lib\site-packages\CADETProcess\optimization\results.py:383, in OptimizationResults.plot_figures(self, show)
    379 self.plot_objectives(
    380     show=show, plot_directory=self.plot_directory
    381 )
    382 if self.optimization_problem.n_variables > 1:
--> 383     self.plot_corner(
    384         show=show, plot_directory=self.plot_directory
    385     )
    387 if self.optimization_problem.n_objectives > 1:
    388     self.plot_pareto(
    389         show=show, plot_directory=self.plot_directory
    390     )

File C:\ProgramData\miniconda3\envs\cadet_env\lib\site-packages\CADETProcess\optimization\results.py:538, in OptimizationResults.plot_corner(self, *args, **kwargs)
    517 """Create a corner plot of the independent variables.
    518 
    519 Parameters
   (...)
    535 
    536 """
    537 try:
--> 538     self.population_all.plot_corner(*args, **kwargs)
    539 except AssertionError:
    540     pass

File C:\ProgramData\miniconda3\envs\cadet_env\lib\site-packages\CADETProcess\optimization\population.py:575, in Population.plot_corner(self, untransformed, show, plot_directory)
    572 x = np.delete(x.transpose(), singular_indices, 0).transpose()
    573 labels = [label for label in labels if label not in singular_labels]
--> 575 fig = corner.corner(
    576     x,
    577     labels=labels,
    578     bins=20,
    579     quantiles=[0.16, 0.5, 0.84],
    580     show_titles=True,
    581     title_kwargs={"fontsize": 20},
    582     title_fmt=".2g",
    583     use_math_text=True,
    584     quiet=True,
    585 )
    586 fig_size = 6*len(labels)
    587 fig.set_size_inches((fig_size, fig_size))

File C:\ProgramData\miniconda3\envs\cadet_env\lib\site-packages\corner\corner.py:280, in corner(data, bins, range, axes_scale, weights, color, hist_bin_factor, smooth, smooth1d, labels, label_kwargs, titles, show_titles, title_quantiles, 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)
    244         logging.warning(
    245             "Please install arviz to use the advanced features of corner"
    246         )
    248     return corner_impl(
    249         data,
    250         bins=bins,
   (...)
    277         **hist2d_kwargs,
    278     )
--> 280 return arviz_corner(
    281     data,
    282     bins=bins,
    283     range=range,
    284     axes_scale=axes_scale,
    285     weights=weights,
    286     color=color,
    287     hist_bin_factor=hist_bin_factor,
    288     smooth=smooth,
    289     smooth1d=smooth1d,
    290     labels=labels,
    291     label_kwargs=label_kwargs,
    292     titles=titles,
    293     show_titles=show_titles,
    294     title_quantiles=title_quantiles,
    295     title_fmt=title_fmt,
    296     title_kwargs=title_kwargs,
    297     truths=truths,
    298     truth_color=truth_color,
    299     scale_hist=scale_hist,
    300     quantiles=quantiles,
    301     verbose=verbose,
    302     fig=fig,
    303     max_n_ticks=max_n_ticks,
    304     top_ticks=top_ticks,
    305     use_math_text=use_math_text,
    306     reverse=reverse,
    307     labelpad=labelpad,
    308     hist_kwargs=hist_kwargs,
    309     group=group,
    310     var_names=var_names,
    311     filter_vars=filter_vars,
    312     coords=coords,
    313     divergences=divergences,
    314     divergences_kwargs=divergences_kwargs,
    315     labeller=labeller,
    316     **hist2d_kwargs,
    317 )

File C:\ProgramData\miniconda3\envs\cadet_env\lib\site-packages\corner\arviz_corner.py:135, in arviz_corner(data, bins, range, axes_scale, 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)
    130     titles = np.concatenate(
    131         [np.asarray(titles[k]).flatten() for k in var_names]
    132     )
    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,
    138     bins=bins,
   (...)
    164     **hist2d_kwargs,
    165 )
    167 # Get diverging draws and combine chains

File <__array_function__ internals>:180, in stack(*args, **kwargs)

File C:\ProgramData\miniconda3\envs\cadet_env\lib\site-packages\numpy\core\shape_base.py:422, 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')
    424 shapes = {arr.shape for arr in arrays}
    425 if len(shapes) != 1:

ValueError: need at least one array to stack
10

Hey Samuel,

Happy to hear! We also enjoyed it a lot! :partying_face:

Regarding your issue, unfortunately, I cannot reproduce it locally. Which version of CADET-Process are you currently running? Note that we recently released v0.8.0 which includes a lot of stability improvements, especially for the optimization module.

To upgrade, run

pip install cadet-process --upgrade

Also, consider adding a transform to your variables. Without one, TrustConstr would quickly be out of bounds for me. I’d recommend log here.

optimization_problem.add_variable(name='flow_sheet.mixer.V', lb=0.01*1e-6, ub=1.5*1e-6, transform='log')
optimization_problem.add_variable(name='flow_sheet.tubing.axial_dispersion', lb=1e-19, ub=1e-10, transform='log')

Finally, consider using another metric, such as Shape. Unfortunately, this also means it’s required to use another Optimizer that can handle multiple objectives, such as U_NSGA3.

Cheers

11

Hey Johannes :v:

I upgraded CADET-Process to v0.8.0 as you recommended, and adding the log-transformation finally led to a successful optimization of the optimization problem.

Thanks for your help!
Samuel

2 Likes