Inverse fit using new CADET

CADET Troubleshooting
1

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;
	
	% Adsorption
	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.quadratic      = zeros(1, model.nComponents);
	model.cubic          = zeros(1, model.nComponents);

	% Sec 1:Loading phase
    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
2

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]
	
	% Adsorption
	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.quadratic      = zeros(2, mLrm.nComponents);
	mIn.cubic          = zeros(2, mLrm.nComponents);

	% Section 1: Loading phase
	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
3

I thank the team for the support.

Regards
Lalita