yemi.py

import numpy as np

import matplotlib.pyplot as plt

from CADETProcess.processModel import ComponentSystem
from CADETProcess.processModel import NoBinding
from CADETProcess.processModel import Source, LumpedRateModelWithPores, Sink
from CADETProcess.processModel import FlowSheet
from CADETProcess.processModel import Process
from CADETProcess.simulator import Cadet
from CADETProcess.comparison import Comparator
from CADETProcess.reference import ReferenceIO
from CADETProcess.optimization import OptimizationProblem
from CADETProcess.optimization import U_NSGA3
from CADETProcess import settings


mydata = np.loadtxt('./75cmph_1.csv', delimiter=',')
time = mydata[:, 0]*60  # Convert time data from minutes to seconds
conductivity = mydata[:, 1]
concentration = ((conductivity/1000)*0.64)/58.44
plt.plot(time, concentration)

mydata = np.loadtxt('./300cmph_1.csv', delimiter=',')
time_300 = mydata[:, 0]*60  # Convert time data from minutes to seconds
conductivity_300 = mydata[:, 1]
concentration_300 = ((conductivity_300/1000)*0.64)/58.44
plt.plot(time_300, concentration_300)

# Set Parameters
inlet_conc = 1000  # 1M NaCl = 1000 mol m-3
feed_flowrate_75 = 0.48  # ml/min
feed_flowrate_300 = 1.92  # ml/min
vol_flowrate_75 = feed_flowrate_75*((1e-6)/60)  # m^3 s-1, Convert feed flow rate units
vol_flowrate_300 = feed_flowrate_300*((1e-6)/60)  # m^3 s-1, Convert feed flow rate units

col_length = 2.5  # cm, column length
col_ID = 0.7  # cm, column inner diameter
col_csa = np.pi*(col_ID/100)**2  # m^2, column cross sectional area
particle_diameter = 5e-5  # m, MSSpCC particle diameter
sim_duration = 840  # s, Equivalent to the load time


#### CADET-Process
component_system = ComponentSystem(1) # One component, i.e. target protein or material
binding_model = NoBinding(component_system, name='nobinding')
feed_unit = Source(component_system, name='feed')
feed_unit.c = [max(concentration)]
column = LumpedRateModelWithPores(component_system, name='column')
column.length = col_length/100  # m
column.diameter = col_ID/100    # m
column.axial_dispersion = 1e-7
column.bed_porosity = 0.4       # -
column.particle_porosity = 0.4  # -
column.particle_radius = particle_diameter/2    # m
column.film_diffusion = [0]     #
column.binding_model = binding_model    # Attach the configured binding model to the column
outlet = Sink(component_system, name='outlet')
flow_sheet = FlowSheet(component_system)
flow_sheet.add_unit(feed_unit)
flow_sheet.add_unit(column)
flow_sheet.add_unit(outlet)
flow_sheet.add_connection(feed_unit, column)
flow_sheet.add_connection(column, outlet)

process_75 = Process(flow_sheet, 'process_75')
process_75.cycle_time = sim_duration
process_75.add_event(
    'feed_on', 'flow_sheet.feed.flow_rate', vol_flowrate_75, 0
)

process_300 = Process(flow_sheet, 'process_300')
process_300.cycle_time = sim_duration
process_300.add_event(
    'feed_on', 'flow_sheet.feed.flow_rate', vol_flowrate_300, 0
)

process_simulator = Cadet()

simulation_results_75 = process_simulator.simulate(process_75)
simulation_results_300 = process_simulator.simulate(process_300)

_ = simulation_results_75.solution.column.outlet.plot()
_ = simulation_results_300.solution.column.outlet.plot()

# %%

# Reference

reference_75 = ReferenceIO('Exp_75', time, concentration)
# reference_75_1 = ReferenceIO('Exp_75_1', time, concentration)

reference_300 = ReferenceIO('Exp_300', time_300, concentration_300)
# reference_300_1 = ReferenceIO('Exp_300_1', time_300, concentration_300)


# Comparator for 75
comparator_75 = Comparator()
comparator_75.add_reference(reference_75)
comparator_75.add_difference_metric('Shape', reference_75, 'column.outlet')

comparator_75.plot_comparison(simulation_results_75)

# Comparator for 300
comparator_300 = Comparator()
comparator_300.add_reference(reference_300)
comparator_300.add_difference_metric('Shape', reference_300, 'column.outlet')

comparator_300.plot_comparison(simulation_results_300)
# %%

# Optimzation Problem
optimization_problem = OptimizationProblem('MSSpCC')
optimization_problem.add_evaluation_object(process_75)
optimization_problem.add_evaluation_object(process_300)
optimization_problem.add_variable(
    name='bed_poros',
    parameter_path='flow_sheet.column.bed_porosity',
    lb=0, ub=1, transform='auto'
)
optimization_problem.add_variable(
    name='axial_dispersion',
    parameter_path='flow_sheet.column.axial_dispersion',
    lb=1e-8, ub=1e-5
)
optimization_problem.add_variable(
    name='particle_porsosity',
    parameter_path='flow_sheet.column.particle_porosity',
    lb=0, ub=1, transform='auto'
)

optimization_problem.add_evaluator(process_simulator)

optimization_problem.add_objective(
    comparator_75,
    n_objectives=comparator_75.n_metrics,
    requires=[process_simulator],
    evaluation_objects=process_75
)
optimization_problem.add_objective(
    comparator_300,
    n_objectives=comparator_300.n_metrics,
    requires=[process_simulator],
    evaluation_objects=process_300
)

settings.working_directory = './yemi'


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


optimization_problem.add_callback(
    plot_75, requires=[process_simulator],
    evaluation_objects=process_75
)


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


optimization_problem.add_callback(
    plot_300, requires=[process_simulator],
    evaluation_objects=process_300
)

optimizer = U_NSGA3()
optimizer.pop_size = 64
optimizer.n_max_gen = 64
optimizer.n_cores = 8


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