# Inverse fit using new CADET

I want to perform the inverse fit for a system consisting of two CSTRs and one column. The sequence of two CSTRs and one column is CSTR-column-CSTR. For the column I want to use lumped rate mode without pores as column model. The inverse fit is performed for three parameters: ka, kd, and column dispersion, in this case. Based on the “parest” folder given in the new CADET I created a MATLAB for performing the inverse fit. I am facing problem in how to incorporate the two CSTRs using “singleCSTR”.

Can I get an example for inverse fit for a system of CSTRs and column.

``````function singleExpSeparateComponents()

% Create simulation and obtain artificial data
sim = createSimulation();
res = sim.run();

% Set model parameters and enable sensitivities
params = cell(3, 1);
params{1} = makeSensitivity([0], {'MCL_KA'},         [1], [-1], [0], [-1]);
params{1} = makeSensitivity([0], {'MCL_KD'},         [2], [-1], [0], [-1]);
params{2} = makeSensitivity([0], {'COL_DISPERSION'}, [3], [-1], [0], [-1]);

sim.setParameters(params, true(3, 1));

% Create fit object
pf = ParameterFit();

% Collect data of first experiment in cell array (one cell per observed component / wavelength)
% Note that time points of the measurements are given in sim.solutionTimes
data = cell(1, 1);
data{1} = res.solution.outlet{1}(:, 2);

% Specify which components are observed in which factor for each observation / wavelength
idxComp = cell(1, 1);
idxComp{1} = [0 1]; % First observation is component 2 (or 1 if read as 0-based) with factor 1
%idxComp{2} = [0 1]; % Second observation is component 3 (or 2 if read as 0-based) with factor 1
%idxComp{3} = [0 0 0 1]; % Third observation is component 4 (or 3 if read as 0-based) with factor 1

% Add the experiment to the fit (unit operation 0 is observed, which is the GRM) including name
% of the experiment and the different observations / wavelengths
pf.addExperiment(data, sim, [0], idxComp, [], 'Frontal run', {'mAb'});

% Specify initial parameters and their lower and upper bounds
initParams = [5, 0.1, 5.75e-8];
loBound = 1e-1 * ones(1, 3);
upBound = 1e2 * ones(1, 3);

% This variable serves as storage for the plot handles returned by the plot function of ParameterFit
% in the OutputFcn below.
plotHd = [];

% Set some options (tolerances, enable Jacobian) and request some output (iteration statistics, plots)
opts = optimset('TolFun', 1e-8, 'TolX', 1e-8, 'MaxIter', 100, 'Jacobian', 'on', 'Diagnostics', 'off', 'Display', 'iter', 'OutputFcn', @progressMonitor);
% Invoke optimizer lsqnonlin
[optimalParams, optimalRes, ~, exitflag] = lsqnonlin(@residual, initParams, loBound, upBound, opts);
success = (exitflag > 0);

function [varargout] = residual(x)
%RESIDUAL Residual function that is passed to lsqnonlin, just forwards to ParameterFit object
varargout = cell(nargout,1);
[varargout{:}] = pf.residualVector(x);
end

function stop = progressMonitor(x, optimValues, state)
%PROGRESSMONITOR Callback function for progress report invoked by lsqnonlin, just calls ParameterFit's plot function
stop = false;
if strcmp(state, 'iter')
% Call plot function and reuse plot handles from previous call (storage outside this function in captured variable plotHd)
plotHd = pf.plot(plotHd, optimValues.iteration, [optimValues.resnorm, optimValues.stepsize]);
end
end

end

function sim = createSimulation()
% Lumped rate model without pores
model = SingleLRM();

model.nComponents = 2; %
model.nCellsColumn = 16; % Attention: This is very low and only used for illustration (short runtime)
model.nBoundStates = ones(model.nComponents, 1);

% Initial conditions
model.initialBulk = [ 176 0.0];
model.initialSolid = [1.2e3 0.0];

% Transport
model.dispersionColumn          = 5.75e-8;
model.interstitialVelocity      = 3e-3;

% Geometry
model.columnLength  = 0.014;
model.porosity      = 0.37;

mLangmuir = LangmuirBinding();
mLangmuir.kineticBinding = true;
mLangmuir.kA         = 2.44;
mLangmuir.kD         = 0.03;
mLangmuir.qMax         = 1.4;
model.bindingModel = mLangmuir;

% Inlet
%mIn = PiecewiseCubicPolyInlet();

model.constant       = zeros(1, model.nComponents);
model.linear         = zeros(1, model.nComponents);
model.cubic          = zeros(1, model.nComponents);

model.constant(1,1)  = 176.0;
model.constant(1,2)  = 0.003013;  % component 1

% Model system
%	mSys = SingleCSTR();
%	mSys.models = [model, mIn];
%	mSys.connectionStartSection = [0];
%	mSys.connections = {[1, 0, -1, -1, 1.0]};

% Configure simulator
sim = Simulator.create();
sim.sectionTimes = [0.0 3776];
sim.sectionContinuity = [];
sim.solutionTimes = linspace(0, 3776, 1000);

% Assign model
sim.model = model;
end
``````

You need to construct a system of unit operations instead of using the `SingleLRM` class.
I’d suggest

1. building the model and running it with some arbitrary (but reasonable) parameters to see if it works;
2. putting the model in a parameter estimation procedure.

In general, I strongly recommend switching to Python. Here, you can use our dedicated parameter fitting framework CADET-Match, which is much more powerful than the Matlab `ParameterFit` class.

This is a fully runnable example for a parameter estimation. The model is set up in `createSimulation()`.

``````function cstrColCstrParEst()

% Create simulation and obtain artificial data
sim = createSimulation();
res = sim.run();

% Set model parameters and enable sensitivities
params = cell(5, 1);
params{1} = makeSensitivity([2], {'MCL_KA'}, [0], [-1], [-1], [0], [-1]);
params{2} = makeSensitivity([2], {'MCL_KA'}, [1], [-1], [-1], [0], [-1]);
params{3} = makeSensitivity([2], {'MCL_KD'}, [0], [-1], [-1], [0], [-1]);
params{4} = makeSensitivity([2], {'MCL_KD'}, [1], [-1], [-1], [0], [-1]);
params{5} = makeSensitivity([2], {'COL_DISPERSION'}, [-1], [-1], [-1], [-1], [-1]);
sim.setParameters(params, true(numel(params), 1));

% Create fit object
pf = ParameterFit();

% Collect data of first experiment in cell array (one cell per observed component / wavelength)
% Note that time points of the measurements are given in sim.solutionTimes
data = cell(2, 1);
data{1} = squeeze(res.solution.outlet{5}(:, 1, 1));
data{2} = squeeze(res.solution.outlet{5}(:, 1, 2));

% Specify which components are observed in which factor for each observation / wavelength
idxComp = cell(2, 1);
idxComp{1} = [1 0]; % First observation is component 1 (or 0 if read as 0-based) with factor 1
idxComp{2} = [0 1]; % Second observation is component 2 (or 1 if read as 0-based) with factor 1

% Add the experiment to the fit (unit operation 4 is observed, which is the system outlet) including name
% of the experiment and the different observations / wavelengths
pf.addExperiment(data, sim, [4 4], idxComp, [], [], [], [], 'Experiment1', {'Comp A', 'Comp B'});

% Specify initial parameters and their lower and upper bounds
initParams = [30, 2, 11, 9, 3e-7];
loBound = [0.1, 0.1, 0.1, 0.1, 1e-8];
upBound = [40, 40, 30, 30, 1e-5];

% This variable serves as storage for the plot handles returned by the plot function of ParameterFit
% in the OutputFcn below.
plotHd = [];

% Set some options (tolerances, enable Jacobian) and request some output (iteration statistics, plots)
opts = optimset('TolFun', 1e-8, 'TolX', 1e-8, 'MaxIter', 100, 'Jacobian', 'on', 'Diagnostics', 'off', 'Display', 'iter', 'OutputFcn', @progressMonitor);
% Invoke optimizer lsqnonlin
[optimalParams, optimalRes, ~, exitflag] = lsqnonlin(@residual, initParams, loBound, upBound, opts);
success = (exitflag > 0);

function [varargout] = residual(x)
%RESIDUAL Residual function that is passed to lsqnonlin, just forwards to ParameterFit object
varargout = cell(nargout,1);
[varargout{:}] = pf.residualVector(x);
end

function stop = progressMonitor(x, optimValues, state)
%PROGRESSMONITOR Callback function for progress report invoked by lsqnonlin, just calls ParameterFit's plot function
stop = false;
if strcmp(state, 'iter')
% Call plot function and reuse plot handles from previous call (storage outside this function in captured variable plotHd)
plotHd = pf.plot(plotHd, optimValues.iteration, [optimValues.resnorm, optimValues.stepsize]);
end
end

end

function sim = createSimulation()
% Step 1: Construct column model
% ===================================================

% General rate model unit operation
mLrm = LumpedRateModelWithoutPores();

% Discretization
mLrm.nComponents = 2;
mLrm.nCellsColumn = 64; % Attention: This is very low and only used for illustration (short runtime)
mLrm.nBoundStates = ones(mLrm.nComponents, 1); % Number of bound states for each component

% Initial conditions (equilibrated empty column)
mLrm.initialBulk = [0.0 0.0]; % [mol / m^3], also used for the particle mobile phase
mLrm.initialSolid = [0.0 0.0]; % [mol / m^3]

% Transport
mLrm.dispersionColumn          = 5.75e-7; % [m^2 / s]
mLrm.interstitialVelocity      = 1; % Forward flow

% Geometry
mLrm.columnLength        = 0.014; % [m]
mLrm.porosity            = 0.8425; % [-]
mLrm.crossSectionArea    = (1.0 / 2)^2 * pi; % [m^2]

mLang = LangmuirBinding();
mLang.kineticBinding = true; % Dynamic binding
mLang.kA         = [35.5 1.59]; % Adsorption rate [m^3 / (mol * s)]
mLang.kD         = [10 10]; % Desorption rate [1 / s]
mLang.qMax       = [4.7 5.29]; % Capacity [mol / m^3]
mLrm.bindingModel = mLang;

% Step 2: Construct inlet unit operation
% ===================================================

% Inlet unit operation
mIn = PiecewiseCubicPolyInlet();
mIn.nComponents = 2;

% Reserve space: nSections x nComponents (a section can be thought of being a
% step in the process, see below)
mIn.constant       = zeros(2, mLrm.nComponents);
mIn.linear         = zeros(2, mLrm.nComponents);
mIn.cubic          = zeros(2, mLrm.nComponents);

mIn.constant(1,1)  = 1.0;   % [mol / m^3] component 3
mIn.constant(1,2)  = 1.0;   % [mol / m^3] component 4

% Step 3: Construct tank unit operations
% ===================================================

% Inlet unit operation
mTank1 = StirredTankModel();
mTank1.nComponents = 2;
mTank1.initialConcentration = [0, 0]; % [mol / m^3]
mTank1.initialVolume = 1e-6; % [m^3]

mTank2 = StirredTankModel();
mTank2.nComponents = 2;
mTank2.initialConcentration = [0, 0]; % [mol / m^3]
mTank2.initialVolume = 1.5e-6; % [m^3]

% Step 4: Construct outlet unit operation
% ===================================================

% Inlet unit operation
mOut = OutletModel();
mOut.nComponents = 2;

% Step 4: Assemble system of unit operations
% ===================================================

% Construct ModelSystem and assign unit operations (order determines IDs)
mSys = ModelSystem();
mSys.models = [mIn, mTank1, mLrm, mTank2, mOut];

flowRate = 1e-4; % [m^3/s]

% Define valve configurations / unit operation connections
% Valve configuration active on entering section 0
mSys.connectionStartSection = [0];
% Connect unit 0 with unit 1 (-1 = all ports, all components), flow rate
% Connect unit 1 with unit 2 (-1 = all ports, all components), flow rate
mSys.connections = {[0, 1, -1, -1, -1, -1, flowRate; ...
1, 2, -1, -1, -1, -1, flowRate; ...
2, 3, -1, -1, -1, -1, flowRate; ...
3, 4, -1, -1, -1, -1, flowRate]};

% Step 5: Create simulator and configure it
% ===================================================

% Construct and configure simulator
sim = Simulator.create();
sim.solutionTimes = linspace(0, 1500, 1001); % [s], time points at which solution is computed

% sectionTimes holds the sections and sectionContinuity indicates whether
% the transition between two adjacent sections is continuous
sim.sectionTimes = [0.0 10.0 1500.0]; % [s]
sim.sectionContinuity = false(1,1);

% Hand model over to simulator
sim.model = mSys;
end
``````

I thank the team for the support.

Regards
Lalita