Figure 9.15:

Model fit to a single adsorption experiment.

Code for Figure 9.15

Text of the GNU GPL.

main.m


  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
% Copyright (C) 2001, James B. Rawlings and John G. Ekerdt
%
% This program is free software; you can redistribute it and/or
% modify it under the terms of the GNU General Public License as
% published by the Free Software Foundation; either version 2, or (at
% your option) any later version.
%
% This program is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
% General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; see the file COPYING.  If not, write to
% the Free Software Foundation, 59 Temple Place - Suite 330, Boston,
% MA 02111-1307, USA.

%
% jbr,  11/28/01
%
%  modified to use paresto, jbr, 6/8/18
%

Hgas1 = [2.659, 27.590, 141.038, ...
	332.613, 546.857, 608.583, 875.353, 991.072, 1080.134, ...
	1142.782, 1203.178, 1256.244]';

Hads1 = [32.31, 40.02, 45.40, 48.49 ...
	50.19, 53.86, 52.77, 53.19, 53.67, 54.15, 54.51, 54.73]';


Hgas2 = [3.755, 36.331, 189.904, ...
	 413.881, 624.665, 817.592]';

Hads2 = [33.17, 40.90, 45.59, 47.72, ...
	 48.74, 49.43]';

Hgas3 = [12.794, 119.331, 339.917, ...
	 545.227, 735.746, 890.642]';

Hads3 = [36.91, 42.90, 45.65, 47.00, ...
	 47.82, 48.60]';

Hgas4 = [7.424, 73.312, 251.177, ...
	 491.691, 707.806, 873.926]';

Hads4 = [39.11, 46.06, 49.87, 51.57, ...
	 52.48, 53.18]';

Hgas5 = [9.176, 72.065, 259.002, ...
	 475.250, 670.034, 818.669]';

Hads5 = [39.52, 46.06, 49.92, 51.86, ...
	 52.96, 53.88]';

Hgas6 = [2.962, 34.350, 176.860, 405.233, ...
	 607.476, 799.723, 946.765]';

Hads6 = [32.70, 40.11, 44.92, 47.40, 48.62, ...
	 49.50, 50.12]';

Hgas7 = [3.159, 38.528, 200.877, 433.924, ...
	 647.630, 820.673]';

Hads7 = [32.42, 39.81, 44.51, 46.79, 48.14, ...
	 48.95]';

Hgas8 = [34.127, 172.316, ...
	 385.633, 616.149, 812.932, 961.191, 1061.967, 1155.829, 1219.299]';

Hads8 = [43.79, 48.97, 51.57, ...
	 53.09, 54.05, 54.67, 55.20, 55.60, 55.93]';

%
% choose a dataset or datasets to fit
%
cadsmeas = [Hads4];
cgmeas = [Hgas4];

Ks0 = 2;
cms0 = 40;
cmp0 = 0;
par0 = [Ks0; cms0; cmp0];

theta0 = [par0; 0];
thetalb = [1e-3; 10; 0; 0];
thetaub = [5; 200; 100; 0];
est_ind = 1:3;

% use a dummy differential state to measure time
algmodel.x = {'time'};
algmodel.z = {'cg', 'cads'};
algmodel.p = {'Ks', 'cms', 'cmp'};
algmodel.d = {'cgmeas', 'cadsmeas'};

algmodel.ode = @(t, y, p) {1};
algmodel.alg = @(t, y, p) {y.cads - p.cmp - ...
			   (p.cms*sqrt(p.Ks)*sqrt(y.cgmeas))/(1 + sqrt(p.Ks)*sqrt(y.cgmeas)), ...
			   y.cg - y.cgmeas};
%algmodel.lsq = @(t, y, p) {(y.cgmeas-y.cg)/10, y.cadsmeas-y.cads};
algmodel.lsq = @(t, y, p) {y.cadsmeas-y.cads};
algmodel.tout = 1:numel(cgmeas);

pe = paresto(algmodel);

est = pe.optimize([cgmeas'; cadsmeas'], theta0, thetalb, thetaub);
theta_opt =  est.theta(est_ind)

theta_conf =  pe.confidence(est, est_ind, 0.95)

nplot = 100;
Xplot = linspace(0,1.05*max(cgmeas),nplot)';
Ks = est.theta(1);
cms = est.theta(2);
cmx = est.theta(3);
Yplot = cms*sqrt(Ks)*sqrt(Xplot)./(1+sqrt(Ks)*sqrt(Xplot)) + cmx;

table1 = [Xplot, Yplot];
table2 = [cgmeas, cadsmeas];
save adsone.dat table1 table2;

if (~ strcmp (getenv ('OMIT_PLOTS'), 'true')) % PLOTTING
plot (table1(:,1), table1(:,2), table2(:,1), table2(:,2), '+');
% TITLE
end % PLOTTING