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 | % 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.
% Revised 7/24/2018
p = struct();
p.deltag1 = -3.72; % kcal/mol
p.deltag2 = -4.49; % kcal/mol
p.T = 400; % K
p.P = 2.5; % atm
p.R = 1.987; % cal/mol/K
p.K1 = exp(-p.deltag1*1000/(p.R*p.T));
p.K2 = exp(-p.deltag2*1000/(p.R*p.T));
% first optimize free energy
% min_z G(z) s.t. x_1,x_2 >0, x_1 + x_2 < 0.5
p.Gval = 0;
p.w0 = [0.2; 0.2];
p.A = [1, 1];
p.b = 0.5;
[w, obj, info] = fmincon(@(w) G(w,p), p.w0, p.A, p.b, [], [], [0;0], [0.5;0.5]);
w
if (info ~= 1)
info, w, obj
end
% next prepare contour plot
Glevels = [0, -1, -2, -2.5, -2.5, -2.53, -2.55 -2.559];
nlev = length(Glevels);
alpha = [0.05, 0.05, 0.075, 0.05, 0.05, 0.05, 0.04, 0.02];
x0vec = [1e-3, 1e-3, 1e-3, .133, .133, .133, .133, .133; ...
.02, .15, .35, .3, .3, .3, .3, .3];
delx = [0.01, 0.01, 0.01, 0.01, -0.01, 0, 0, 0; ...
0 0, 0, 0, 0, 0, 0, 0];
sfinal = [9, 9, 9, 9, 9, 9, 9, 5];
ztable(1:nlev) = {NaN(2,1)};
for i = 1:nlev
p.Gval = Glevels(i);
% find a point on the contour corresponding to s=0
ext0 = x0vec(:,i);
unk = [2];
x0 = ext0(unk);
[x, fval, info] = fsolve(@(x) Gcontour(x,p,ext0,unk), x0);
fsolve_failed = info <= 0;
infovec(i) = info;
fcnvec(:,i) = Gcontour(x,p,ext0,unk);
ext = ext0;
ext(unk) = x;
zlast = ext;
extinit(:,i) = ext;
% continuation in arclength
dels = alpha(i);
svec = dels:dels:sfinal(i);
narc =length(svec);
z0 = ext + delx(:,i);
unk = [1 2];
for j = 1:narc
s = svec(j);
[z, fval, info] = fsolve(@(z) inflate(z,zlast,dels,p,ext0,unk),z0);
fsolve_ok = info > 0;
% stop when contour completed
if (z(1) > 0 && z(2) > 0 && 0.5-(z(1)+z(2)) > 0 && fsolve_ok)
ztable{i}(:,j) = z;
% project initial guess along the arc for next point on contour
z0 = 2*z - zlast;
zlast = z;
else
break
end
end
ztable(i) = {[ext, ztable{i}]'};
end
feasbound = [0, 0.5; 0.5, 0]; % store boundary of feasible region
save contour_pentane.dat ztable feasbound;
%plot results
if (~ strcmp (getenv ('OMIT_PLOTS'), 'true')) % PLOTTING
hold on;
plot([0.5,0], [0, 0.5]); % one leg of feasible region
for i = 1:nlev
plot(ztable{i}(:,1), ztable{i}(:,2));
end
% TITLE
end % PLOTTING
|
