* authors: Moll Glass and Alexander Mitsos, RWTH Aachen University, AVT - Aachener Verfahrenstechnik, Process Systems Engineering, Aachen, Germany * 07.09.2017 * You may use this code as long as you acknowledge the source. * the purpose of this file is to regress parameters of the proposed Gibbs free energy model of Ruszczynski et al (J. Chem. Eng. Data, 2017) * as opposed to the authors' regression method, we use a global solver and simultaneously regress all parameters * this allows for the best possible fit upon convergence of the respective solver * imposing the constraint c > - 2*T enforces partial convexity of the submodels * OctOH/H2O - example of a broad miscibility gap * BARON converges to an excellent fit, c1 and c2 are at upper bounds of stability stated by Ruszczynski et al (J. Chem. Eng. Data, 2017) * data source: Dallos, A.; Liszi, J. J. Chem. Thermodyn., 1995, 27, 447-8 (Liquid + liquid) equilibria of (octan-1-ol + water) at temperatures from 2 88.15 K to 323.15 K * note that Baker's criterion is met, and the correct number of phases can be guaranteed, but this is not the case for the number of phase splits $ontext calculations were run on a 64 bit windows 7 (service pack 1) system with an Intel Core i3-3225 3.30 GHz processor GAMS 24.8.5, using BARON 17.4.1 with default options the reported results are: modelstatus 2.00 solverstatus 1.00 fuBestPointFound 0.000004143226315 fuFinalLowerBound -0.000005856773685 a1 0.815888381918774 a2 5.830584534646562 b1 141.998625850303100 b2 1067.881389907787000 c1 -0.000010000000000 c2 -0.000010000000000 x_phase1 0.729838447249950 x_phase2 0.000052667715932 x_phase1 0.727559016895547 x_phase2 0.000055929076139 x_phase1 0.725337817402873 x_phase2 0.000059270784793 x_phase1 0.723172718928993 x_phase2 0.000062689983347 x_phase1 0.721061689619338 x_phase2 0.000066183805010 x_phase1 0.719002790413587 x_phase2 0.000069749384951 x_phase1 0.716994170140257 x_phase2 0.000073383869414 x_phase1 0.715034060886258 x_phase2 0.000077084423821 CPU time 0.040000000000000 measurement data No. Liquid1 mole fraction OCTANOL Temperature (K) Pressure (N/sqm) Liquid2 mole fraction OCTANOL 1 0,729 288,15 101325 5,4E-05 2 0,727 293,15 101325 5,6E-05 3 0,726 298,15 101325 5,9E-05 4 0,724 303,15 101325 6,2E-05 5 0,722 308,15 101325 6,5E-05 6 0,719 313,15 101325 6,9E-05 7 0,717 318,15 101325 7,4E-05 8 0,714 323,15 101325 7,8E-05 $offtext set experiments /exp1*exp8/; * mole fractions in mol/mol variable x_phase1(experiments); x_phase1.LO(experiments)=1e-6; x_phase1.L(experiments)=0.5; x_phase1.UP(experiments)=1-1e-6; variable x_phase2(experiments); x_phase2.LO(experiments)=1e-6; x_phase2.L(experiments)=0.5; x_phase2.UP(experiments)=1-1e-6; * parameters variable a1; a1.LO=-10; a1.UP=10; a1.L=0; variable a2; a2.LO=-10; a2.UP=10; a2.L=0; variable b1; b1.LO=-10; b1.L=0; b1.UP=4000; variable b2; b2.LO=-10; b2.L=0; b2.UP=4000; * the lower bounds satisfy -2T (-576 K) to ensure partial convexity * the upper bounds strictly satisfy condition stated by Ruszczynski et al (J. Chem. Eng. Data, 2017) variable c1; c1.LO=-500; c1.L=0; c1.UP=-1e-5; variable c2; c2.LO=-500; c2.L=0; c2.UP=-1e-5; * temperature in K parameter T(experiments); T('exp1')= 288.15; T('exp2')= 293.15; T('exp3')= 298.15; T('exp4')= 303.15; T('exp5')= 308.15; T('exp6')= 313.15; T('exp7')= 318.15; T('exp8')= 323.15; * isopotential (Equation (15)) in Ruszczynski et al (J. Chem. Eng. Data, 2017) * since we have constraints on c-parameters, Baker is satisfied if isopotential is met Equation eq_isopotential_species1(experiments); eq_isopotential_species1(experiments) .. log(x_phase1(experiments))-c1/T(experiments)*sqr(1-x_phase1(experiments)) =E= log(x_phase2(experiments))+c2/T(experiments)*(2*x_phase2(experiments)-sqr(x_phase2(experiments)))+ a2 + b2/T(experiments) ; Equation eq_isopotential_species2(experiments); eq_isopotential_species2(experiments) .. log(1-x_phase1(experiments)) + c1/T(experiments)*(1-sqr(x_phase1(experiments))) + a1 + b1/T(experiments) =E= log(1-x_phase2(experiments))- c2/T(experiments) * sqr(x_phase2(experiments)) ; * we use the same objective function (Equation 20) as Ruszczynski et al (J. Chem. Eng. Data, 2017) * However, we recommend to use different weights for the individual phases since the miscbility gap is very broad variable objective; Equation eq_objective; eq_objective .. objective =E= 0 +sqr(0.729-x_phase1('exp1')) +10000*sqr(5.4e-5-x_phase2('exp1')) +sqr(0.727-x_phase1('exp2')) +10000*sqr(5.6e-5-x_phase2('exp2')) +sqr(0.726-x_phase1('exp3')) +10000*sqr(5.9e-5-x_phase2('exp3')) +sqr(0.724-x_phase1('exp4')) +10000*sqr(6.2e-5-x_phase2('exp4')) +sqr(0.722-x_phase1('exp5')) +10000*sqr(6.5e-5-x_phase2('exp5')) +sqr(0.719-x_phase1('exp6')) +10000*sqr(6.9e-5-x_phase2('exp6')) +sqr(0.717-x_phase1('exp7')) +10000*sqr(7.4e-5-x_phase2('exp7')) +sqr(0.714-x_phase1('exp8')) +10000*sqr(7.8e-5-x_phase2('exp8')) ; model mymodel /all/; options optca=0.00001; options optcr=0; options nlp=baron; solve mymodel using nlp minimizing objective; file output /regression_res.txt/; put output; put 'modelstatus ',mymodel.modelstat/; put 'solverstatus ',mymodel.SolveStat/; put 'fuBestPointFound ',mymodel.objVal:25:15/ ; put 'fuFinalLowerBound ',mymodel.objEst:25:15/ ; put 'a1',a1.L:25:15/; put 'a2',a2.L:25:15/; put 'b1',b1.L:25:15/; put 'b2',b2.L:25:15/; put 'c1',c1.L:25:15/; put 'c2',c2.L:25:15/; loop(experiments, put 'x_phase1', x_phase1.L(experiments):25:15/; put 'x_phase2', x_phase2.L(experiments):25:15/; put /; ); put 'CPU time', mymodel.resUsd:25:15/;