Stoichiometric coefficient estimation in CADET

CADET Troubleshooting
1

Hello!

I am working on a project involving a reaction model with the goal of estimating stoichiometric coefficients to fit experimental data of a chemical reaction. My primary objective is to adjust the reaction model’s parameters so that simulated results align closely with experimental observations.

The reaction I want to adjust is the following one:
glucose + glycerol + NH3 → product

Experimental data is available for the concentration of these components over time and I am using the MassActionLaw class. To estimate the stoichiometric coefficients I am formulating the problem as an optimization problem, but I am having trouble correctly implementing this. Is it possible to solve this kind of problem using CADET? So far, I have tried implementing the following code, but without success:

import pandas as pd
from CADETProcess.processModel import ComponentSystem, MassActionLaw, Cstr, FlowSheet, Process
from CADETProcess.simulator import Cadet
from CADETProcess.reference import ReferenceIO
from CADETProcess.comparison import Comparator
from CADETProcess.optimization import OptimizationProblem, U_NSGA3

data_path = "experimental_data.xlsx"
data = pd.read_excel(data_path, index_col=0)

component_system = ComponentSystem(["glucose", "glycerol", "NH3", "CDW"])

reaction_system = MassActionLaw(component_system)
reaction_system.add_reaction(
    indices=[0, 1, 2, 3],
    coefficients=[-1.0, -1.0, -1.0, 1.0],
    k_fwd=1.0, 
    k_bwd=0.0
)

reactor = Cstr(component_system, "reactor")
reactor.bulk_reaction_model = reaction_system
reactor.V = 1.0
reactor.c = [1.0, 1.0, 1.0, 0.0]

flow_sheet = FlowSheet(component_system, "flow_sheet")
flow_sheet.add_unit(reactor)

process = Process(flow_sheet, "reaction_fermentation")
process.cycle_time = 180

reference = ReferenceIO("experimental data", data.index, data, component_system=component_system)

comparator = Comparator("fermentation_comparator")
comparator.add_reference(reference)
comparator.add_difference_metric("RMSE", reference, "reactor.outlet", components=["glucose", "glycerol", "NH3", "CDW"])

optimization_problem = OptimizationProblem("fermentation_optimization")
optimization_problem.add_evaluation_object(process)

optimization_problem.add_variable(
    name="stoichiometric_coefficient",
    parameter_path="flow_sheet.reactor.bulk_reaction_model.stoich",
    indices=(0),
    lb=-5.0,
    ub=5.0
)

simulator = Cadet()
optimization_problem.add_evaluator(simulator)

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

optimizer = U_NSGA3()
optimizer.n_cores = 4
optimizer.pop_size = 32 
optimizer.n_max_gen = 20

optimization_results = optimizer.optimize(optimization_problem)

print(optimization_results)

optimization_results.plot_convergence()
optimization_results.plot_objectives()

Thank you very much in advance.

Kind regards,
Ruoxian

2

Hello Ruoxian,

that should absolutely be possible but there appears to be a bug in CADET-Process (or a misconfiguration in the optimizationProblem that I’m not spotting). I’ve opened a bug report on GitHub and will update you once we figured it out.

Best wishes,
Ron

3

Hi,

I believe I have fixed the bug. Please try installing the corresponding branch with the following command

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

and report back whether this solves your problem.

1 Like
4

Thanks to the quick fix by @j.schmoelder this should work now.

A few closing remarks:

  1. Thank you for bringing this to our attention and for sharing code that can reproduce the error. There are, as always, things that can make it even easier for developers to track down bugs:
    a. shared code should ideally be an MRE (minimal reproducible example) without requiring non-included data. In your code you use an additional xlsx file data_path = "experimental_data.xlsx" that we don’t have access to. In this case, removing the references to the file was easy enough because the bug was reproducible without starting the optimizer, but for future bug-reports please ensure the code runs without additional input.
    b. Including the verbatim error output or a detailed description of where/ how it fails is very helpful for developers.
  2. The stochiometry matrix is a 2D matrix of (#components, #reactions). That becomes more evident if you add a second reaction and then look at the output of reaction_system.stoich. To be explicit, it is therefore advisable to specify the full index when setting up the optimization variable like so:

# optimize stochiometric coefficient of the 2nd component (index 1) 
#  of the first reaction (index 0):
optimization_problem.add_variable(
    name="stoichiometric_coefficient",
    parameter_path="flow_sheet.reactor.bulk_reaction_model.stoich",
    indices=(1, 0),  # index is (#component, #reaction)
    lb=-5.0,
    ub=5.0
)

(In this example I’ve chosen to optimize the stochiometric coefficient of component 2 (Glycerol) to ensure different indexes between component and reaction.) Hope that makes sense.

If you have further questions, please don’t hesitate to reach out. We also have the Office Hours where you can drop by to chat with us.

1 Like
5

Hi,

this is now merged into dev. Hopefully, it fixes your problem. However, we would still appreciate your feedback.

Best regards

Jo