diff --git a/SMB.m b/SMB.m index ebb3c56..3aeda6e 100644 --- a/SMB.m +++ b/SMB.m @@ -41,12 +41,43 @@ model = ModelGRM(); model.nComponents = opt.nComponents; - model.kineticBindingModel = false; - model.bindingModel = LinearBinding(); % if you want to change the equilibrium isotherm -% Adsorption parameters - model.bindingParameters.LIN_KA = opt.KA; - model.bindingParameters.LIN_KD = opt.KD; +% if you want to change the equilibrium isotherm + if strcmp(opt.BindingModel, 'LinearBinding') + + model.kineticBindingModel = false; + model.bindingModel = LinearBinding(); + +% Adsorption parameters + model.bindingParameters.LIN_KA = opt.KA; + model.bindingParameters.LIN_KD = opt.KD; + + elseif strcmp(opt.BindingModel, 'MultiComponentLangmuirBinding') + + model.kineticBindingModel = true; + model.bindingModel = MultiComponentLangmuirBinding(); + + model.bindingParameters.MCL_KA = opt.KA; + model.bindingParameters.MCL_KD = opt.KD; + model.bindingParameters.MCL_QMAX = opt.QMAX; + + elseif strcmp(opt.BindingModel, 'MultiComponentBiLangmuirBinding') + + model.kineticBindingModel = true; + model.bindingModel = MultiComponentBiLangmuirBinding(); + + model.bindingParameters.MCL_KA1 = opt.KA(1); + model.bindingParameters.MCL_KD1 = opt.KD(1); + model.bindingParameters.MCL_QMAX1 = opt.QMAX(1); + model.bindingParameters.MCL_KA2 = opt.KA(2); + model.bindingParameters.MCL_KD2 = opt.KD(2); + model.bindingParameters.MCL_QMAX2 = opt.QMAX(2); + + elseif strcmp(opt.BindingModel, 'StericMassAction') + + error('%s: it is not available yet.', opt.BindingModel); + + end if nargin >= 3 && ~isempty(lastState) model.initialState = lastState; @@ -156,7 +187,7 @@ % | | | | | % Zone II(b) | | Zone II(d/c) | % | | | | | -% Ext1 —— Zone IV(d) | Ext1 —— Zone IV(g/h) +% Ext1 -- Zone IV(d) | Ext1 -- Zone IV(g/h) % | | | | | % Zone I(a) | | Zone I(b/a) | % | | | | | @@ -171,7 +202,7 @@ % | | | | | % Zone II(f/e/d) | | Zone II(h/g/f/e) | % | | | | | -% Ext1 —— Zone IV(j/k/l) | Ext1 —— Zone IV(m/n/o/p) +% Ext1 -- Zone IV(j/k/l) | Ext1 -- Zone IV(m/n/o/p) % | | | | | % Zone I(c/b/a) | | Zone I(d/c/b/a) | % | | | | | @@ -1431,7 +1462,7 @@ end -% Please be quite careful, which components is used for statistics (change them with comp_ext_ID or comp_raf_ID) +% Please be quite careful, which component is used for statistics (change them with comp_ext_ID or comp_raf_ID) if opt.nComponents == 2 % Extract ports Purity_extract = trapz(plotData{1,position_ext}.outlet.time, plotData{1,position_ext}.outlet.concentration(:,opt.comp_ext_ID)) /... @@ -1490,7 +1521,7 @@ position_ext1 = 18; position_ext2 = 14; position_raf = 6; end -% Please be quite careful, which components is used for statistics (change them with comp_ext_ID or comp_raf_ID) +% Please be quite careful, which component is used for statistics (change them with comp_ext_ID or comp_raf_ID) if opt.nComponents == 3 % Extract ports Purity_extract1 = trapz(plotData{1,position_ext1}.outlet.time, plotData{1,position_ext1}.outlet.concentration(:,opt.comp_ext1_ID)) /... @@ -2161,7 +2192,7 @@ function plotFigures(opt, plotData) end - end % if opt.nZone == 4/ opt.nZone == 5 + end % if opt.nZone == 4 / opt.nZone == 5 end % if opt.enableDebug @@ -2181,4 +2212,4 @@ function plotFigures(opt, plotData) % Institute: Forschungszentrum Juelich GmbH, IBG-1, Juelich, Germany. % % All rights reserved. Please see the license of CADET. -% ============================================================================= \ No newline at end of file +% ============================================================================= diff --git a/SMBOptimization.m b/SMBOptimization.m index d3766a9..29c38c4 100644 --- a/SMBOptimization.m +++ b/SMBOptimization.m @@ -41,29 +41,29 @@ function SMBOptimization() % The set of the parameters which are optimized params = struct('columnLength',[], 'switch',[], 'recycle',[], 'feed',[], 'desorbent',[], 'extract1',[], 'extract2',[]); - + % Select one method and make it true (correspondingly the rest methods false) - optimization_method.Particle_Swarm_Optimization = true; optimization_method.Differential_Evolution = false; - optimization_method.Metropolis_Adjusted_Differential_Evolution = false; + optimization_method.Particle_Swarm_Optimization = true; optimization_method.Deterministic_algorithm_fmincon = false; + optimization_method.Metropolis_Adjusted_Differential_Evolution = false; if isfield(optimization_method, 'Particle_Swarm_Optimization') ... && optimization_method.Particle_Swarm_Optimization - + OptAlgorithms.Particle_Swarm_Optimization(params); - + elseif isfield(optimization_method, 'Differential_Evolution') ... && optimization_method.Differential_Evolution - + OptAlgorithms.Differential_Evolution(params); - + elseif isfield(optimization_method, 'Metropolis_Adjusted_Differential_Evolution') ... && optimization_method.Metropolis_Adjusted_Differential_Evolution - + OptAlgorithms.Metropolis_Adjusted_Differential_Evolution(params); - + % elseif isfield(optimization_method, 'Riemann_Manifold_Metropolis_Adjusted_Langevin') ... % && optimization_method.Riemann_Manifold_Metropolis_Adjusted_Langevin % @@ -78,10 +78,10 @@ function SMBOptimization() loBound = [0.20, 150, 8.0e-7, 0.9e-7, 0.7e-7, 1.0e-7]; upBound = [0.30, 230, 10e-7, 2.0e-7, 2.0e-7, 2.0e-7]; - + options = optimoptions('fmincon', 'Algorithm', 'interior-point', 'Display', 'iter',... 'TolX',1e-6,'TolCon',1e-6,'TolFun',1e-6,'MaxIter',500); - + try [SMBparams, fval, exitflag, output, ~, grad] = fmincon( @simulatedMovingBed, ... initParams, [],[],[],[], loBound, upBound, [], options); @@ -89,16 +89,16 @@ function SMBOptimization() disp('Errors in the MATLAB build-in optimizer: fmincon. \n Please check your input parameters and run again. \n'); disp('The message from fmincon: %s \n', exception.message); end - + fprintf('Minimum: %g, Parameters:[%g| %g| %g| %g| %g| %g] \n', fval, SMBparams); - + else - + warning('The method you selected is not provided in this programme'); - + end - - + + end % ============================================================================= % SMB - The Simulated Moving Bed Chromatography for separation of diff --git a/examples/Forward/getParameters_binary_ModiCon.m b/examples/Forward/getParameters_binary_ModiCon.m index 1d7ea53..707f6af 100644 --- a/examples/Forward/getParameters_binary_ModiCon.m +++ b/examples/Forward/getParameters_binary_ModiCon.m @@ -35,11 +35,9 @@ opt.enableDebug = true; % set it true when you want to see the figures opt.nZone = 4; % 4-zone for binary separation, 5-zone for ternary separation opt.nColumn = 4; % 4,8,12,16 -column cases are available -% opt.nColumn = 8; -% opt.nColumn = 12; -% opt.nColumn = 16; % Binding: Linear Binding isotherm + opt.BindingModel = 'LinearBinding'; opt.nComponents = 2; opt.KA = [5.72 7.7]; % [comp_A, comp_B], A for raffinate, B for extract opt.KD = [1, 1]; diff --git a/examples/Forward/getParameters_binary_case1.m b/examples/Forward/getParameters_binary_case1.m index 779bb76..e991c65 100644 --- a/examples/Forward/getParameters_binary_case1.m +++ b/examples/Forward/getParameters_binary_case1.m @@ -35,11 +35,9 @@ opt.enableDebug = true; opt.nZone = 4; % 4-zone for binary separation, 5-zone for ternary separation opt.nColumn = 4; % 4,8,12,16- column cases are available -% opt.nColumn = 8; -% opt.nColumn = 12; -% opt.nColumn = 16; % Binding: Linear Binding isotherm + opt.BindingModel = 'LinearBinding'; opt.nComponents = 2; opt.KA = [5.72 7.7]; % [comp_A, comp_B], A for raffinate, B for extract opt.KD = [1, 1]; diff --git a/examples/Forward/getParameters_binary_case2.m b/examples/Forward/getParameters_binary_case2.m index 6ac0913..1faf468 100644 --- a/examples/Forward/getParameters_binary_case2.m +++ b/examples/Forward/getParameters_binary_case2.m @@ -35,10 +35,9 @@ opt.enableDebug = true; opt.nZone = 4; % 4-zone for binary separation, 5-zone for ternary separation opt.nColumn = 8; % 4,8,12,16- column cases are available -% opt.nColumn = 12; -% opt.nColumn = 16; % Binding: Linear Binding isotherm + opt.BindingModel = 'LinearBinding'; opt.nComponents = 2; opt.KA = [0.28 0.54]; % [comp_A, comp_B], A for raffinate, B for extract opt.KD = [1, 1]; diff --git a/examples/Forward/getParameters_binary_ternary.m b/examples/Forward/getParameters_binary_ternary.m index 02d5e77..173f3b9 100644 --- a/examples/Forward/getParameters_binary_ternary.m +++ b/examples/Forward/getParameters_binary_ternary.m @@ -35,11 +35,9 @@ opt.enableDebug = true; opt.nZone = 4; % 4-zone for binary separation, 5-zone for ternary separation opt.nColumn = 8; % 4,8,12,16- column cases are available -% opt.nColumn = 8; -% opt.nColumn = 12; -% opt.nColumn = 16; % Binding: Linear Binding isotherm + opt.BindingModel = 'LinearBinding'; opt.nComponents = 3; opt.KA = [0.23, 0.28, 0.61]; % [comp_A comp_B comp_C], A for raffinate, B C for extract opt.KD = [1, 1, 1]; @@ -89,7 +87,6 @@ Feed.concentration(1:end,i) = (concentrationFeed(i) / opt.molMass(i)); end - end % ============================================================================= % SMB - The Simulated Moving Bed Chromatography for separation of diff --git a/examples/Forward/getParameters_ternary_case1.m b/examples/Forward/getParameters_ternary_case1.m index e6c037f..f275fd4 100644 --- a/examples/Forward/getParameters_ternary_case1.m +++ b/examples/Forward/getParameters_ternary_case1.m @@ -35,10 +35,10 @@ opt.enableDebug = true; % set it true when you want to see the figures opt.nZone = 5; % 5-zone for ternary separation - opt.nColumn = 5; % 5-column case is available so far -% opt.nColumn = 10; + opt.nColumn = 10; % 5, 10 -column case is available so far % Binding: Linear Binding isotherm + opt.BindingModel = 'LinearBinding'; opt.nComponents = 3; opt.KA = [3.15, 7.4, 23]; % [comp_A, comp_B, comp_C], A,B for raffinate, C for extract opt.KD = [1, 1, 1]; % K_A < K_B < K_C diff --git a/examples/Forward/getParameters_ternary_case2.m b/examples/Forward/getParameters_ternary_case2.m index 002db1b..a07beba 100644 --- a/examples/Forward/getParameters_ternary_case2.m +++ b/examples/Forward/getParameters_ternary_case2.m @@ -35,10 +35,10 @@ opt.enableDebug = true; % set it true when you want to see the figures opt.nZone = 5; % 5-zone for ternary separation -% opt.nColumn = 5; % 5-column case is available so far - opt.nColumn = 10; + opt.nColumn = 10; % 5, 10 -column case is available % Binding: Linear Binding isotherm + opt.BindingModel = 'LinearBinding'; opt.nComponents = 3; opt.KA = [5.34, 6.80, 11.20]; % [comp_A, comp_B, comp_C], A,B for raffinate, C for extract opt.KD = [1, 1, 1]; % K_A < K_B < K_C diff --git a/examples/Optimization/getParameters_binary.m b/examples/Optimization/getParameters_binary.m index b84a76c..8d41636 100644 --- a/examples/Optimization/getParameters_binary.m +++ b/examples/Optimization/getParameters_binary.m @@ -17,7 +17,7 @@ valueAssign = struct('columnLength',ParSwarm(1), 'switch',ParSwarm(2), 'recycle',ParSwarm(3),... 'feed',ParSwarm(4), 'desorbent',ParSwarm(5), 'extract',ParSwarm(6)); - + % The parameter setting for simulator opt.tolIter = 1e-3; % tolerance of the SMB stopping criterion opt.nMaxIter = 1000; % the maximum iteration step in SMB @@ -27,28 +27,26 @@ opt.ABSTOL = 1e-9; % tolerance of CADET stopping criterion opt.INIT_STEP_SIZE = 1e-14; % refer your to CADET manual opt.MAX_STEPS = 5e6; % the maximum iteration step in CADET - + % The parameter setting for the SMB opt.switch = valueAssign.switch; % switching time opt.timePoints = 1000; % the observed time-points opt.Purity_extract_limit = 0.99; % used for constructing constraints opt.Purity_raffinate_limit = 0.99; % used for constructing constraints opt.Penalty_factor = 10; % penalty factor in penalty function - + opt.enableDebug = false; % set it false if you are using the optimizer opt.nZone = 4; opt.nColumn = 4; % 4,8,12,16 -column cases are available -% opt.nColumn = 8; -% opt.nColumn = 12; -% opt.nColumn = 16; % Binding: Linear Binding isotherm + opt.BindingModel = 'LinearBinding'; opt.nComponents = 2; opt.KA = [5.72 7.7]; % [comp_A, comp_B], A for raffinate, B for extract opt.KD = [1, 1]; opt.comp_raf_ID = 1; % the target component withdrawn from the raffinate ports opt.comp_ext_ID = 2; % the target component withdrawn from the extract ports - + % Transport opt.dispersionColumn = 3.8148e-20; % D_{ax} opt.filmDiffusion = [100 100]; % K_{eff} @@ -79,11 +77,11 @@ interstVelocity.raffinate = flowRate.raffinate / (crossArea*opt.porosityColumn); % m/s interstVelocity.desorbent = flowRate.desorbent / (crossArea*opt.porosityColumn); % m/s interstVelocity.extract = flowRate.extract / (crossArea*opt.porosityColumn); % m/s - + concentrationFeed = [0.55, 0.55]; % g/m^3 [concentration_compA, concentration_compB] opt.molMass = [180.16, 180.16]; % The molar mass of each components opt.yLim = max(concentrationFeed ./ opt.molMass); % the magnitude for plotting - + % Feed concentration setup Feed.time = linspace(0, opt.switch, opt.timePoints); Feed.concentration = zeros(length(Feed.time), opt.nComponents); @@ -91,7 +89,7 @@ for i = 1:opt.nComponents Feed.concentration(1:end,i) = (concentrationFeed(i) / opt.molMass(i)); end - + end % ============================================================================= % SMB - The Simulated Moving Bed Chromatography for separation of diff --git a/examples/Optimization/getParameters_ternary.m b/examples/Optimization/getParameters_ternary.m index 2fc4b05..68a6438 100644 --- a/examples/Optimization/getParameters_ternary.m +++ b/examples/Optimization/getParameters_ternary.m @@ -17,7 +17,7 @@ valueAssign = struct('columnLength',ParSwarm(1), 'switch',ParSwarm(2), 'recycle',ParSwarm(3),... 'feed',ParSwarm(4), 'desorbent',ParSwarm(5), 'extract',ParSwarm(6)); - + % The parameter setting for simulator opt.tolIter = 1e-3; % tolerance of the SMB stopping criterion opt.nMaxIter = 1000; % the maximum iteration step in SMB @@ -27,7 +27,7 @@ opt.ABSTOL = 1e-9; % tolerance of CADET stopping criterion opt.INIT_STEP_SIZE = 1e-14; % refer your to CADET manual opt.MAX_STEPS = 5e6; % the maximum iteration step in CADET - + % The parameter setting for the SMB opt.switch = valueAssign.switch; % s % switching time opt.timePoints = 1000; % the observed time-points @@ -35,20 +35,20 @@ opt.Purity_extract2_limit = 0.69; % used for constructing constraints opt.Purity_raffinate_limit = 0.99; % used for constructing constraints opt.Penalty_factor = 10; % penalty factor in penalty function - + opt.enableDebug = false; % set it false if you are using the optimizer opt.nZone = 5; opt.nColumn = 5; % 4,8,12,16 -column cases are available -% opt.nColumn = 10; % Binding: Linear Binding isotherm + opt.BindingModel = 'LinearBinding'; opt.nComponents = 3; opt.KA = [3.15, 7.4, 23]; % [comp_A, comp_B], A for raffinate, B for extract opt.KD = [1, 1, 1]; opt.comp_raf_ID = 1; % the target component withdrawn from the raffinate ports opt.comp_ext1_ID = 3; % the target component withdrawn from the extract ports opt.comp_ext2_ID = 2; % the target component withdrawn from the extract ports - + % Transport opt.dispersionColumn = 3.8148e-10; % D_{ax} opt.filmDiffusion = [5.0e-5, 2.5e-5, 5.0e-5]; % K_f @@ -82,11 +82,11 @@ interstVelocity.desorbent = flowRate.desorbent / (crossArea*opt.porosityColumn); % m/s interstVelocity.extract1 = flowRate.extract1 / (crossArea*opt.porosityColumn); % m/s interstVelocity.extract2 = flowRate.extract2 / (crossArea*opt.porosityColumn); % m/s - + concentrationFeed = [1.0, 1.0, 1.0]; % g/m^3 [concentration_compA, concentration_compB] opt.molMass = [227.217, 267.24, 251.24192]; % The molar mass of each components opt.yLim = max(concentrationFeed ./ opt.molMass); % the magnitude for plotting - + % Feed concentration setup Feed.time = linspace(0, opt.switch, opt.timePoints); Feed.concentration = zeros(length(Feed.time), opt.nComponents); @@ -94,7 +94,7 @@ for i = 1:opt.nComponents Feed.concentration(1:end,i) = (concentrationFeed(i) / opt.molMass(i)); end - + end % ============================================================================= % SMB - The Simulated Moving Bed Chromatography for separation of diff --git a/getParameters.m b/getParameters.m new file mode 100644 index 0000000..1faf468 --- /dev/null +++ b/getParameters.m @@ -0,0 +1,100 @@ +function [opt, interstVelocity, Feed] = getParameters(varargin) +% Case 2, a eight-column demonstration case + +% ============================================================================= +% This is the function to input all the necessary data for simulation +% +% Returns: +% 1. opt stands for options, which involves the parameter settings +% for the algorithm, the binding isotherm, and the model equations +% +% 2. interstVelocity is calculated from flowrate of each column and inlet. +% interstitial_velocity = flow_rate / (across_area * porosity_Column) +% +% 3. Feed initializes the injection concentration +% ============================================================================= + + +% The parameter setting for simulator + opt.tolIter = 1e-4; + opt.nMaxIter = 1000; + opt.nThreads = 8; + opt.nCellsColumn = 40; + opt.nCellsParticle = 1; + opt.ABSTOL = 1e-10; + opt.INIT_STEP_SIZE = 1e-14; + opt.MAX_STEPS = 5e6; + +% The parameter settting for the SMB + opt.switch = 1552; + opt.timePoints = 1000; + opt.Purity_extract_limit = 0.99; + opt.Purity_raffinate_limit = 0.99; + opt.Penalty_factor = 10; + + opt.enableDebug = true; + opt.nZone = 4; % 4-zone for binary separation, 5-zone for ternary separation + opt.nColumn = 8; % 4,8,12,16- column cases are available + +% Binding: Linear Binding isotherm + opt.BindingModel = 'LinearBinding'; + opt.nComponents = 2; + opt.KA = [0.28 0.54]; % [comp_A, comp_B], A for raffinate, B for extract + opt.KD = [1, 1]; + opt.comp_raf_ID = 1; % the target component withdrawn from the raffinate ports + opt.comp_ext_ID = 2; % the target component withdrawn from the extract ports + +% Transport + opt.dispersionColumn = 3.8148e-6; % + opt.filmDiffusion = [5e-5 5e-5]; % unknown + opt.diffusionParticle = [1.6e4 1.6e4]; % unknown + opt.diffusionParticleSurface = [0.0 0.0]; + +% Geometry + opt.columnLength = 53.6e-2; % m + opt.columnDiameter = 2.60e-2; % m + opt.particleRadius = 0.325e-2 /2; % m + opt.porosityColumn = 0.38; + opt.porosityParticle = 0.00001; % unknown + +% Parameter units transformation +% The flow rate of Zone I was defined as the recycle flow rate + crossArea = pi * (opt.columnDiameter/2)^2; % m^2 + flowRate.recycle = 0.1395e-6; % m^3/s + flowRate.feed = 0.02e-6 ; % m^3/s + flowRate.raffinate = 0.0266e-6; % m^3/s + flowRate.desorbent = 0.0414e-6; % m^3/s + flowRate.extract = 0.0348e-6; % m^3/s + opt.flowRate_extract = flowRate.extract; + opt.flowRate_raffinate = flowRate.raffinate; + +% Interstitial velocity = flow_rate / (across_area * opt.porosityColumn) + interstVelocity.recycle = flowRate.recycle / (crossArea*opt.porosityColumn); % m/s + interstVelocity.feed = flowRate.feed / (crossArea*opt.porosityColumn); % m/s + interstVelocity.raffinate = flowRate.raffinate / (crossArea*opt.porosityColumn); % m/s + interstVelocity.desorbent = flowRate.desorbent / (crossArea*opt.porosityColumn); % m/s + interstVelocity.extract = flowRate.extract / (crossArea*opt.porosityColumn); % m/s + + concentrationFeed = [0.5, 0.5]; % g/m^3 [concentration_compA, concentration_compB] + opt.molMass = [180.16, 180.16]; + opt.yLim = max(concentrationFeed ./ opt.molMass); + +% Feed concentration setup + Feed.time = linspace(0, opt.switch, opt.timePoints); + Feed.concentration = zeros(length(Feed.time), opt.nComponents); + + for i = 1:opt.nComponents + Feed.concentration(1:end,i) = (concentrationFeed(i) / opt.molMass(i)); + end + +end +% ============================================================================= +% SMB - The Simulated Moving Bed Chromatography for separation of +% target compounds, either binary or ternary. +% +% Author: QiaoLe He E-mail: q.he@fz-juelich.de +% +% Institute: Forschungszentrum Juelich GmbH, IBG-1, Juelich, Germany. +% +% All rights reserved. Please see the license of CADET. +% ============================================================================= \ No newline at end of file diff --git a/isSMBupdateAvailable.m b/isSMBupdateAvailable.m index df4af48..5ec4c05 100644 --- a/isSMBupdateAvailable.m +++ b/isSMBupdateAvailable.m @@ -6,18 +6,20 @@ localPath = fileparts(mfilename('fullpath')); + fprintf('Adding %s to MATLAB PATH\n', localPath); + path(localPath, path); fileID = fopen([localPath filesep 'version.txt']); - + currentVersion = []; try currentVersion = fgets(fileID); end - + if isempty(currentVersion) fprintf('The local version.txt file has been deleted. \n'); return; end - + splitted = textscan(currentVersion, '%s', 'delimiter', '.'); currentVersion = splitted{1}; @@ -25,12 +27,12 @@ try stableVersion = urlread('https://raw.githubusercontent.com/modsim/CADET-SMB/master/version.txt'); end - + if isempty(stableVersion) fprintf('The internet connection is not available now. \n'); return; end - + stableVersionString = stableVersion; splitted = textscan(stableVersion, '%s', 'delimiter', '.'); stableVersion = splitted{1}; @@ -42,7 +44,7 @@ break; end end - + fprintf('The newest version is now installed, Version %s', stableVersionString); end diff --git a/version.txt b/version.txt index 5154b3f..95e3ba8 100644 --- a/version.txt +++ b/version.txt @@ -1 +1 @@ -2.6 +2.5