combined different Binding Models; prompt to the version 2.6

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.
-% =============================================================================
+% =============================================================================

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,27 +78,27 @@ 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);
         catch exception
             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

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];

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];

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];

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

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

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

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,19 +77,19 @@
     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);
 
     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

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,28 +27,28 @@
     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
     opt.Purity_extract1_limit  = 0.95;  % used for constructing constraints
     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,19 +82,19 @@
     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);
 
     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