diff --git a/.idea/dictionaries/Titus.xml b/.idea/dictionaries/Titus.xml index 54115dc..7ca660c 100644 --- a/.idea/dictionaries/Titus.xml +++ b/.idea/dictionaries/Titus.xml @@ -7,6 +7,7 @@ equilibrate extractant kmol + lmse ndarray pred reeps diff --git a/.idea/workspace.xml b/.idea/workspace.xml index 5ba0cb0..ee250e0 100644 --- a/.idea/workspace.xml +++ b/.idea/workspace.xml @@ -2,11 +2,11 @@ - - - + + + - - + + - + @@ -208,7 +208,14 @@ @@ -227,7 +234,8 @@ - @@ -246,10 +254,10 @@ - + - + @@ -290,10 +298,10 @@ - + - + diff --git a/data/xmls/twophase.xml b/data/xmls/twophase.xml index cb9e012..28f9edb 100644 --- a/data/xmls/twophase.xml +++ b/data/xmls/twophase.xml @@ -50,7 +50,7 @@ 298.14999999999998 - -4848809.807683906 + -4704703.645715787 1117.965 0.0 diff --git a/reeps.py b/reeps.py index f28314c..727fd23 100644 --- a/reeps.py +++ b/reeps.py @@ -1,4 +1,4 @@ -import time +from datetime import datetime import cantera as ct import pandas as pd import numpy as np @@ -7,7 +7,8 @@ from scipy.optimize import minimize import xml.etree.ElementTree as ET import seaborn as sns import matplotlib.pyplot as plt - +import shutil +import copy sns.set() @@ -17,8 +18,7 @@ class REEPS: Returns parameters for GEM :param exp_csv_filename: (str) csv file name with experimental data :param phases_xml_filename: (str) xml file with parameters for equilibrium calc - :param x_guess: (float) guess for multiplier variable - :param h_guess: (float) initial guess standard enthalpy (J/kmol) + :param opt_dict: (dict) optimize info {species:{thermo_prop:guess} :param phase_names: (list) names of phases in xml file :param aq_solvent_name: (str) name of aqueous solvent in xml file :param extractant_name: (str) name of extractant in xml file @@ -34,8 +34,7 @@ class REEPS: def __init__(self, exp_csv_filename, phases_xml_filename, - x_guess, - h_guess, + opt_dict, phase_names, aq_solvent_name, extractant_name, @@ -48,8 +47,7 @@ class REEPS: ): self._exp_csv_filename = exp_csv_filename self._phases_xml_filename = phases_xml_filename - self._x_guess = x_guess - self._h_guess = h_guess + self._opt_dict = opt_dict self._phase_names = phase_names self._aq_solvent_name = aq_solvent_name self._extractant_name = extractant_name @@ -60,6 +58,8 @@ class REEPS: self._extractant_rho = extractant_rho self._diluant_rho = diluant_rho + self._temp_xml_filename = "temp.xml" + shutil.copyfile(phases_xml_filename, self._temp_xml_filename) self._phases = ct.import_phases(phases_xml_filename, phase_names) self._exp_df = pd.read_csv(self._exp_csv_filename) @@ -90,18 +90,11 @@ class REEPS: self.set_in_moles(feed_vol=1) return None - def get_x_guess(self) -> float: - return self._x_guess - - def set_x_guess(self, x_guess): - self._x_guess = x_guess - return None - - def get_h_guess(self) -> float: - return self._h_guess + def get_opt_dict(self) -> dict: + return self._opt_dict - def set_h_guess(self, h_guess): - self._h_guess = h_guess + def set_opt_dict(self, opt_dict): + self._opt_dict = opt_dict return None def get_aq_solvent_name(self) -> str: @@ -161,7 +154,8 @@ class REEPS: return None def set_in_moles(self, feed_vol): - """:param feed_vol: (float) feed volume of mixture (g/L)""" + """Function that initializes mole fractions + :param feed_vol: (float) feed volume of mixture (g/L)""" phases_copy = self._phases.copy() exp_df = self._exp_df.copy() solvent_name = self._aq_solvent_name @@ -244,8 +238,10 @@ class REEPS: return self._in_moles def objective(self, x): - h = x[0] * self._h_guess - phases_copy = self._phases.copy() + """Log mean squared error between measured and predicted data + :param x: (list) thermo properties varied to minimize LMSE""" + temp_xml_filename = self._temp_xml_filename + phase_names = self._phase_names aq_ind = self._aq_ind org_ind = self._org_ind complex_name = self._complex_name @@ -253,18 +249,16 @@ class REEPS: in_moles = self._in_moles exp_df = self._exp_df - complex_index = phases_copy[org_ind].species_index(complex_name) - complex_species = phases_copy[org_ind].species(complex_index) - - new_thermo_coeffs = complex_species.thermo.coeffs - new_thermo_coeffs[1] = h - complex_species.thermo = ct.ConstantCp( - complex_species.thermo.min_temp, - complex_species.thermo.max_temp, - complex_species.thermo.reference_pressure, - new_thermo_coeffs) - phases_copy[org_ind].modify_species(complex_index, complex_species) + x = np.array(x) + opt_dict = copy.deepcopy(self._opt_dict) + i = 0 + for species_name in opt_dict.keys(): + for thermo_prop in opt_dict[species_name].keys(): + opt_dict[species_name][thermo_prop] *= x[i] + i += 1 + self.update_xml(opt_dict, temp_xml_filename) + phases_copy = ct.import_phases(temp_xml_filename, phase_names) mix = ct.Mixture(phases_copy) pred = [] for row in in_moles.values: @@ -285,30 +279,48 @@ class REEPS: Returns estimated complex enthalpy in J/mol :param kwargs: (dict) parameters and options for scipy.minimize """ - x_guess = self._x_guess - h_guess = self._h_guess + opt_dict = copy.deepcopy(self._opt_dict) + # x_guess = [] + i = 0 + for species_name in opt_dict.keys(): + for _ in opt_dict[species_name].keys(): + # x_guess.append(opt_dict[species_name][thermo_prop]) + i += 1 + x_guess = np.ones(i) res = minimize(self.objective, x_guess, **kwargs) - pred_complex_enthalpy = res.x[0] * h_guess / 1000 # J/mol - return pred_complex_enthalpy + i = 0 + for species_name in opt_dict.keys(): + for thermo_prop in opt_dict[species_name].keys(): + opt_dict[species_name][thermo_prop] *= res.x[i] + i += 1 + return opt_dict def update_xml(self, - complex_enthalpy, - phases_xml_filename=None, - complex_name=None): - """updates xml file with new complex_enthalpy""" + info_dict, + phases_xml_filename=None): + """updates xml file with info_dict + :param info_dict: (dict) info in {species_names:{thermo_prop:val}} + :param phases_xml_filename: (str) xml filename if editing other xml + """ if phases_xml_filename is None: phases_xml_filename = self._phases_xml_filename - if complex_name is None: - complex_name = self._complex_name tree = ET.parse(phases_xml_filename) root = tree.getroot() # Update xml file - for species in root.iter('species'): - if species.attrib['name'] is complex_name: - for h0 in species.iter('h0'): - h0.text = str(complex_enthalpy) - h0.set('updated', 'Updated at ' + time.strftime('%X')) + for species_name in info_dict.keys(): + for thermo_prop in info_dict[species_name].keys(): + for species in root.iter('species'): + if species.attrib['name'] == species_name: + for changed_prop in species.iter(thermo_prop): + changed_prop.text = str( + info_dict[species_name][thermo_prop]) + now = datetime.now() + changed_prop.set('updated', + 'Updated at {0}:{1} {2}-{3}-{4}' + .format(now.hour, now.minute, + now.month, now.day, + now.year)) tree.write(phases_xml_filename) diff --git a/tests/one_comp_settings.txt b/tests/one_comp_settings.txt deleted file mode 100644 index 4d2182a..0000000 --- a/tests/one_comp_settings.txt +++ /dev/null @@ -1 +0,0 @@ -{"exp_csv_filename": "../data/csvs/exp_data.csv", "phases_xml_filename": "../data/xmls/twophase.xml", "x_guess": 0.96, "h_guess": -4856609000.0, "phase_names": ["HCl_electrolyte", "PC88A_liquid"], "aq_solvent_name": "H2O(L)", "extractant_name": "(HA)2(org)", "diluant_name": "dodecane", "complex_name": "Nd(H(A)2)3(org)", "rare_earth_ion_name": "Nd+++", "aq_solvent_rho": 1000.0, "extractant_rho": 960.0, "diluant_rho": 750.0} \ No newline at end of file diff --git a/tests/test_reeps.py b/tests/test_reeps.py index bc672f0..911f1e6 100644 --- a/tests/test_reeps.py +++ b/tests/test_reeps.py @@ -5,8 +5,9 @@ sys.path.append('../') from reeps import REEPS -with open('one_comp_settings.txt') as file: +with open('one_ree_settings.txt') as file: testing_params = json.load(file) + beaker = REEPS(**testing_params) minimizer_kwargs = {"method": 'SLSQP', @@ -15,5 +16,7 @@ minimizer_kwargs = {"method": 'SLSQP', "options": {'disp': True, 'maxiter': 1000, 'ftol': 1e-6}} est_enthalpy = beaker.fit(minimizer_kwargs) print(est_enthalpy) +# info_dict = {"Nd(H(A)2)3(org)": {"h0": est_enthalpy}} +# beaker.update_xml(est_enthalpy) beaker.parity_plot()