Enquiry about add_parameter_sensitivity

CADET Troubleshooting
1

Hello all,

I am trying to use the sensitivity analysis feature of cadet-process. Specifically, in my code, I would like to obtain the partial derivative of the concentration at the outlet of tube1 c_out(L,t) with respect to the concentration at the inlet of tube1 c_in(0,t), i.e., the Jacobian matrix ( s = ∂c_out(L,t)/∂c_in(0,t) ). I would like to ask if this can be achieved using the add_parameter_sensitivity function? Additionally, is my current code configuration correct? It seems that the Jacobian matrix Iobtained is all zeros.

Thanks for any help and advice!

import numpy as np
import matplotlib.pyplot as plt
from CADETProcess.simulator import Cadet
from CADETProcess.processModel import ComponentSystem
from CADETProcess.processModel import StericMassAction
from CADETProcess.processModel import Inlet
from CADETProcess.processModel import Cstr
from CADETProcess.processModel import LumpedRateModelWithoutPores, LRMDiscretizationFV
from CADETProcess.processModel import Outlet
from CADETProcess.processModel import FlowSheet
from CADETProcess.processModel import Process
import pandas as pd

def data_generation(CV, c_pro):
    component_system = ComponentSystem()
    component_system.add_component('salt')
    component_system.add_component('Protein1')

    simulator = Cadet()
    lrm_dis_fv = LRMDiscretizationFV()
    lrm_dis_fv.use_analytic_jacobian = False

    Tubediameter = 0.0005
    cross_section_area = np.pi * (0.025 / 2) ** 2

    Capacity = 560  # mM
    kde = 100
    keq = 0.489
    v = 5.94
    sigma = 50
    kad = kde * keq

    Buffer_A = 1  # mM
    Buffer_B = 1001  # mM

    # Feed
    feed = Inlet(component_system, 'feed')
    feed.c = [Buffer_A, c_pro]  # mM

    # Buffer
    buffer = Inlet(component_system, 'buffer')
    buffer.c = [Buffer_A, 0]

    # Tube1
    tube1 = LumpedRateModelWithoutPores(component_system, 'tube1')
    tube1Volume = 10e-8
    tube1.diameter = Tubediameter
    tube1.length = tube1Volume / (np.pi * (Tubediameter / 2) ** 2)
    tube1.total_porosity = 1
    tube1.axial_dispersion = 1e-9  # m^2/s
    tube1.c = [Buffer_A, 0]
    tube1.discretization.use_analytic_jacobian = False

    column1 = LumpedRateModelWithoutPores(component_system, 'column1')
    column1.length = 0.002
    column1.cross_section_area = cross_section_area
    column1.diameter = np.sqrt(4 * column1.cross_section_area / np.pi)
    column1.total_porosity = 0.724
    column1.axial_dispersion = 1e-7
    column1.discretization.use_analytic_jacobian = False

    binding_model = StericMassAction(component_system, 'SMA')
    binding_model.is_kinetic = True
    binding_model.adsorption_rate = [0, kad]
    binding_model.desorption_rate = [0, kde]
    binding_model.characteristic_charge = [0, v]
    binding_model.steric_factor = [0, sigma]
    binding_model.capacity = Capacity
    column1.binding_model = binding_model
    column1.c = [Buffer_A, 0]
    column1.q = [Capacity, 0]

    cstr1 = Cstr(component_system, 'cstr1')
    cstr1.V = 10e-8
    cstr1.c = [Buffer_A, 0]

    tube2 = LumpedRateModelWithoutPores(component_system, 'tube2')
    tube2Volume = 30e-8
    tube2.diameter = Tubediameter
    tube2.length = tube2Volume / (np.pi * (Tubediameter / 2) ** 2)
    tube2.total_porosity = 1
    tube2.axial_dispersion = 1e-9  # m^2/s
    tube2.c = [Buffer_A, 0]
    tube2.discretization.use_analytic_jacobian = False

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

    # Flow Sheet
    flow_sheet = FlowSheet(component_system, 'flow_sheet')
    flow_sheet.add_unit(feed)
    flow_sheet.add_unit(buffer)
    flow_sheet.add_unit(tube1)
    flow_sheet.add_unit(cstr1)
    flow_sheet.add_unit(tube2)
    flow_sheet.add_unit(column1)
    flow_sheet.add_unit(outlet)

    flow_sheet.add_connection(feed, tube1)
    flow_sheet.add_connection(buffer, tube1)
    flow_sheet.add_connection(tube1, column1)


    flow_sheet.add_connection(column1, cstr1)
    flow_sheet.add_connection(cstr1, tube2)
    flow_sheet.add_connection(tube2, outlet)
    process = Process(flow_sheet, 'process')

    flow_rate = 8.6e-6/60

    injection_time = 5e-6 / flow_rate
    Wash_time = 60

    process.add_event(
        'feed_on',
        'flow_sheet.feed.flow_rate',
        flow_rate,
        time=0
    )

    process.add_event(
        'Wash_Start',
        'flow_sheet.feed.flow_rate',
        0,
        injection_time
    )

    process.add_event(
        'Elution',
        'flow_sheet.buffer.flow_rate',
        flow_rate
    )

    process.add_event_dependency('Elution', ['Wash_Start'])

    process.add_event(
        'load', 'flow_sheet.buffer.c',
        [Buffer_A, 0], 0)

    c_wash = np.array([Buffer_A, 0])
    c_elute = np.array([Buffer_B, 0])
    gradient_slope = (c_elute - c_wash) / (CV*6)
    c_gradient_poly = np.array(list(zip(c_wash, gradient_slope)))

    process.add_event(
        'elute_start',
        'flow_sheet.buffer.c',
        c_gradient_poly, injection_time + Wash_time
    )


    process.cycle_time = injection_time + Wash_time + CV*6  # s


    process.add_parameter_sensitivity(
        'feed.c',
        name='feed.c',
        components=['Protein1'],
        polynomial_coefficients='0',
        section_indices=0,
    )

    results = simulator.simulate(process)
    sen = results.sensitivity['feed.c'].tube1.outlet.solution

    return results



c_pro = 0.11

results = data_generation(20, c_pro)
t_SIM = results.solution.outlet.outlet.time
c_salt_SIM = results.solution.outlet.outlet.solution[:, 0]
c_pro_SIM = results.solution.outlet.outlet.solution[:, 1]

fig, ax1 = plt.subplots(figsize=(10, 6))
ax1.set_xlabel('time s')
ax1.set_ylabel('Protein / mM', color='tab:blue')
ax1.plot(t_SIM, c_pro_SIM, '-', label='Pro Sim.', color='C0', linewidth=3)

ax2 = ax1.twinx()
ax2.set_ylabel('Salt / mM', color='tab:green')
ax2.plot(t_SIM, c_salt_SIM, '-', color='black', label='Gradient Sim.', linewidth=3)
ax1.set_xlim(0, 240)
ax1.set_ylim(0, 0.3)
plt.tight_layout()
plt.show()
2

Hey Yang,
thank you for posting this!
We’re short-staffed this week, but I’ll take a look at some point.

Is this the smallest possible process for which this is happening? If not, it would be helpful to identify an even smaller process with only one column and one component that interests you.
Also, could you please let me know which version of CADET-Process you are using?

Best wishes,
Antonia

3

Hi Antonia,

Thank you very much for your help.
I identified a smaller process that only includes feed, buffer, columu, and outlet, but the Jacobian matrix is still all zero. I am using CADET-Process V0.11.1 and CADET V5.0.3.

Best regards,
Yang

4

Thank you so much! Would it be possible to share the smaller example as well?
That would make it easier for me to understand the issue.

5

Yes. Thank you for your help!

import numpy as np
import matplotlib.pyplot as plt
from CADETProcess.simulator import Cadet
from CADETProcess.processModel import ComponentSystem
from CADETProcess.processModel import StericMassAction
from CADETProcess.processModel import Inlet
from CADETProcess.processModel import Cstr
from CADETProcess.processModel import LumpedRateModelWithoutPores, LRMDiscretizationFV
from CADETProcess.processModel import Outlet
from CADETProcess.processModel import FlowSheet
from CADETProcess.processModel import Process
import pandas as pd

def data_generation(CV, c_pro):
    component_system = ComponentSystem()
    component_system.add_component('salt')
    component_system.add_component('Protein1')

    simulator = Cadet()
    lrm_dis_fv = LRMDiscretizationFV()
    lrm_dis_fv.use_analytic_jacobian = False

    Tubediameter = 0.0005
    cross_section_area = np.pi * (0.025 / 2) ** 2

    Capacity = 560  # mM
    kde = 100
    keq = 0.489
    v = 5.94
    sigma = 50
    kad = kde * keq

    Buffer_A = 1  # mM
    Buffer_B = 1001  # mM

    # Feed
    feed = Inlet(component_system, 'feed')
    feed.c = [Buffer_A, c_pro]  # mM

    # Buffer
    buffer = Inlet(component_system, 'buffer')
    buffer.c = [Buffer_A, 0]

    column1 = LumpedRateModelWithoutPores(component_system, 'column1')
    column1.length = 0.002
    column1.cross_section_area = cross_section_area
    column1.diameter = np.sqrt(4 * column1.cross_section_area / np.pi)
    column1.total_porosity = 0.724
    column1.axial_dispersion = 1e-7
    column1.discretization.use_analytic_jacobian = False

    binding_model = StericMassAction(component_system, 'SMA')
    binding_model.is_kinetic = True
    binding_model.adsorption_rate = [0, kad]
    binding_model.desorption_rate = [0, kde]
    binding_model.characteristic_charge = [0, v]
    binding_model.steric_factor = [0, sigma]
    binding_model.capacity = Capacity
    column1.binding_model = binding_model
    column1.c = [Buffer_A, 0]
    column1.q = [Capacity, 0]

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

    # Flow Sheet
    flow_sheet = FlowSheet(component_system, 'flow_sheet')
    flow_sheet.add_unit(feed)
    flow_sheet.add_unit(buffer)
    flow_sheet.add_unit(column1)
    flow_sheet.add_unit(outlet)

    flow_sheet.add_connection(feed, column1)
    flow_sheet.add_connection(buffer, column1)
    flow_sheet.add_connection(column1, outlet)

    process = Process(flow_sheet, 'process')

    flow_rate = 8.6e-6/60

    injection_time = 5e-6 / flow_rate
    Wash_time = 60

    process.add_event(
        'feed_on',
        'flow_sheet.feed.flow_rate',
        flow_rate,
        time=0
    )

    process.add_event(
        'Wash_Start',
        'flow_sheet.feed.flow_rate',
        0,
        injection_time
    )

    process.add_event(
        'Elution',
        'flow_sheet.buffer.flow_rate',
        flow_rate
    )

    process.add_event_dependency('Elution', ['Wash_Start'])

    process.add_event(
        'load', 'flow_sheet.buffer.c',
        [Buffer_A, 0], 0)

    c_wash = np.array([Buffer_A, 0])
    c_elute = np.array([Buffer_B, 0])
    gradient_slope = (c_elute - c_wash) / (CV*6)
    c_gradient_poly = np.array(list(zip(c_wash, gradient_slope)))

    process.add_event(
        'elute_start',
        'flow_sheet.buffer.c',
        c_gradient_poly, injection_time + Wash_time
    )


    process.cycle_time = injection_time + Wash_time + CV*6  # s


    process.add_parameter_sensitivity(
        'feed.c',
        name='feed.c',
        components=['Protein1'],
        polynomial_coefficients='0',
        section_indices=0,
    )

    results = simulator.simulate(process)
    sen = results.sensitivity['feed.c'].column1.outlet.solution

    return results



c_pro = 0.11

results = data_generation(20, c_pro)
t_SIM = results.solution.outlet.outlet.time
c_salt_SIM = results.solution.outlet.outlet.solution[:, 0]
c_pro_SIM = results.solution.outlet.outlet.solution[:, 1]

fig, ax1 = plt.subplots(figsize=(10, 6))
ax1.set_xlabel('time s')
ax1.set_ylabel('Protein / mM', color='tab:blue')
ax1.plot(t_SIM, c_pro_SIM, '-', label='Pro Sim.', color='C0', linewidth=3)

ax2 = ax1.twinx()
ax2.set_ylabel('Salt / mM', color='tab:green')
ax2.plot(t_SIM, c_salt_SIM, '-', color='black', label='Gradient Sim.', linewidth=3)
ax1.set_xlim(0, 240)
ax1.set_ylim(0, 0.3)
plt.tight_layout()
plt.show()
1 Like
6

Hey Yang,

Thank you!
I have some good news and some not-as-good news for you.

First, the good news: Your script looks good, except for two minor issues. I got it to work with these changes:

  1. In the ‘Elution’ event, you need to specify the time at which the event should occur.

  2. You can delete the “load” event:

process.add_event(
“load”, “flow_sheet.buffer.c”,
        [Buffer_A, 0], 0)

because these are your initial conditions.

There is the config, which works for me.

Now for the not-as-good news: I haven’t been able to get the sensitivities to work yet; there also seems to be something wrong with the name of the parameter.

I will get back to you on this.

Best regards
Antonia

1 Like
7

Hi Antonia,

Thank you for pointing out the issue! I have corrected it.

Regards,
Yang

8

This topic was automatically closed 2 days after the last reply. New replies are no longer allowed.