diff --git a/InstallOnMac.sh b/InstallOnMac.sh old mode 100644 new mode 100755 index a719529..def1729 --- a/InstallOnMac.sh +++ b/InstallOnMac.sh @@ -5,9 +5,13 @@ echo -------------------------------------------------- echo Enter environment name: read env_name echo ------------------------------ -echo Creating pandas environment... +echo Creating conda environment... echo ------------------------------ conda create -n $env_name python=3.9 --file requirements.txt -c conda-forge +echo ------------------------------ +echo Activating conda environment... +echo ------------------------------ +conda init conda activate $env_name echo ---------------------- echo Installing PharmaPy... diff --git a/PharmaPy/Kinetics.py b/PharmaPy/Kinetics.py index a39b0e9..3fa4a98 100644 --- a/PharmaPy/Kinetics.py +++ b/PharmaPy/Kinetics.py @@ -239,6 +239,8 @@ def __init__(self, path, k_params, ea_params, rxn_list=None, stoich_matrix = get_stoich(di, partic_species) else: stoich_matrix = np.atleast_2d(stoich_matrix) + if partic_species is None: + raise PharmaPyTypeError('Please provide a participating species list when using a stoichiometric matrix.') perm_idx = get_permutation_indexes(name_species, partic_species) stoich_matrix = stoich_matrix[:, perm_idx] diff --git a/PharmaPy/ParamEstim.py b/PharmaPy/ParamEstim.py index 9f1a285..c6bc5e6 100644 --- a/PharmaPy/ParamEstim.py +++ b/PharmaPy/ParamEstim.py @@ -22,7 +22,11 @@ from PharmaPy.LevMarq import levenberg_marquardt from PharmaPy.Commons import plot_sens, reorder_sens -from cyipopt import minimize_ipopt +try: + from cyipopt import minimize_ipopt + have_cyipopt = True +except ImportError: + have_cyipopt = False linestyles = cycle(['-', '--', '-.', ':']) @@ -615,6 +619,8 @@ def optimize_fn(self, optim_options=None, simulate=False, verbose=True, **optim_options) elif method == 'IPOPT': + if not have_cyipopt: + raise ImportError('cyipopt is an optional import. Please install cyipopt to use IPOPT as a solver for parameter estimation. conda install -c conda-forge cyipopt') if optim_options is None: optim_options = {'print_level': int(verbose) * 5} else: diff --git a/PharmaPy/Reactors.py b/PharmaPy/Reactors.py index 66e38c0..3a079fa 100644 --- a/PharmaPy/Reactors.py +++ b/PharmaPy/Reactors.py @@ -19,7 +19,7 @@ from PharmaPy.CheckModule import check_modeling_objects import numpy as np -from numpy.core.umath_tests import inner1d +# from numpy.core.umath_tests import inner1d import matplotlib.pyplot as plt from matplotlib.ticker import AutoMinorLocator from matplotlib.animation import FuncAnimation @@ -125,10 +125,10 @@ class _BaseReactor: f(time) = state_value return_sens : bool (optional, default = True) whether or not the paramest_wrapper method should return - the sensitivity system along with the concentratio profiles. + the sensitivity system along with the concentration profiles. Use False if you want the parameter estimation platform to estimate the sensitivity system using finite differences - state_events : lsit of dict(s) + state_events : list of dict(s) list of dictionaries, each one containing the specification of a state event """ @@ -506,7 +506,7 @@ def plot_profiles(self, pick_comp=None, **fig_kwargs): can be either the name of a species (str) or the index of the species (int). The default is None. **fig_kwargs : keyword arguments to plt.subplots() - named arguments passed to the plotting functions. A yypical field + named arguments passed to the plotting functions. A typical field is 'figsize', passed as a (width, height) tuple. Returns @@ -605,7 +605,7 @@ class BatchReactor(_BaseReactor): 'coil', 'bath'] return_sens : bool (optional, default = True) whether or not the paramest_wrapper method should return - the sensitivity system along with the concentratio profiles. + the sensitivity system along with the concentration profiles. Use False if you want the parameter estimation platform to estimate the sensitivity system using finite differences """ @@ -688,7 +688,10 @@ def energy_balances(self, time, mole_conc, vol, temp, temp_ht, inputs, delta_hrxn=deltah_rxn) # Balance terms (W) - source_term = -inner1d(deltah_rxn, rates) * vol*1000 # vol in L + #source_term = -inner1d(deltah_rxn, rates) * vol * 1000 # vol in L + # TODO: Check if this is correct + # source_term = -np.inner(deltah_rxn, rates) * vol * 1000 # vol in L + source_term = -(deltah_rxn * rates).sum(axis=1) * vol * 1000 # vol in L if heat_prof: if 'temp' in self.controls.keys(): @@ -942,7 +945,7 @@ class CSTR(_BaseReactor): 'coil', 'bath'] return_sens : bool (optional, default = True) whether or not the paramest_wrapper method should return - the sensitivity system along with the concentratio profiles. + the sensitivity system along with the concentration profiles. Use False if you want the parameter estimation platform to estimate the sensitivity system using finite differences """ @@ -1045,7 +1048,10 @@ def energy_balances(self, time, mole_conc, vol, temp, temp_ht, inputs, flow_term = inlet_flow * (h_in - h_temp) # W # Balance terms (W) - convert vol to L - source_term = -inner1d(deltah_rxn, rates) * vol * 1000 + # source_term = -inner1d(deltah_rxn, rates) * vol * 1000 + # TODO: Check if this is correct + # source_term = -np.dot(deltah_rxn, rates) * vol * 1000 # vol in L + source_term = -(deltah_rxn * rates).sum(axis=1) * vol * 1000 # vol in L if heat_prof: if self.isothermal: @@ -1242,7 +1248,7 @@ class SemibatchReactor(CSTR): 'coil', 'bath'] return_sens : bool (optional, default = True) whether or not the paramest_wrapper method should return - the sensitivity system along with the concentratio profiles. + the sensitivity system along with the concentration profiles. Use False if you want the parameter estimation platform to estimate the sensitivity system using finite differences """ @@ -1442,7 +1448,7 @@ class PlugFlowReactor(_BaseReactor): 'coil', 'bath'] return_sens : bool (optional, default = True) whether or not the paramest_wrapper method should return - the sensitivity system along with the concentratio profiles. + the sensitivity system along with the concentration profiles. Use False if you want the parameter estimation platform to estimate the sensitivity system using finite differences """ @@ -1550,7 +1556,10 @@ def energy_steady(self, conc, temp): delta_hrxn=deltah_rxn) # ---------- Balance terms (W) - source_term = -inner1d(deltah_rxn, rates) * 1000 # W/m**3 + # source_term = -inner1d(deltah_rxn, rates) * 1000 # W/m**3 + # TODO: Check if this is correct + # source_term = -np.dot(deltah_rxn, rates) * 1000 # W / m**3 + source_term = -(deltah_rxn * rates).sum(axis=1) * 1000 # W / m**3 if self.adiabatic: heat_transfer = 0 @@ -1645,7 +1654,11 @@ def energy_balances(self, time, mole_conc, vol_diff, temp, flow_in, rate_i, _, cp_j = self.Liquid_1.getCpPure(temp) # Volumetric heat capacity - cp_vol = inner1d(cp_j, mole_conc) * 1000 # J/m**3/K + # cp_vol = inner1d(cp_j, mole_conc) * 1000 # J/m**3/K + # TODO: Check if this is correct + # cp_vol = -np.dot(cp_j, mole_conc) * 1000 # J/m**3/K + cp_vol = (cp_j * mole_conc).sum(axis=1) * 1000 # vol in L + # Heat of reaction delta_href = self.Kinetics.delta_hrxn @@ -1656,7 +1669,10 @@ def energy_balances(self, time, mole_conc, vol_diff, temp, flow_in, rate_i, stoich, temp, self.mask_species, delta_href, tref_hrxn) # J/mol # ---------- Balance terms (W) - source_term = -inner1d(deltah_rxn, rate_i * 1000) # W/m**3 + # source_term = -inner1d(deltah_rxn, rate_i * 1000) # W/m**3 + # TODO: Check if this is correct + # source_term = -np.dot(deltah_rxn, rate_i * 1000) # W/m**3 + source_term = -(deltah_rxn * rate_i * 1000).sum(axis=1) # W/m**3 temp_diff = np.diff(temp) flow_term = -flow_in * temp_diff / vol_diff # K/s @@ -1910,7 +1926,7 @@ def retrieve_results(self, time, states): def plot_steady(self, fig_size=None, title=None): fig, axes = plt.subplots(1, 2, figsize=fig_size) - # Concentration + # Concentrationn axes[0].plot(self.volPosition, self.concProfSteady) axes[0].set_ylabel('$C_j$ (mol/L)') axes[0].legend(self.name_species) diff --git a/PharmaPy/Results.py b/PharmaPy/Results.py index affecdf..d3b772a 100644 --- a/PharmaPy/Results.py +++ b/PharmaPy/Results.py @@ -103,11 +103,11 @@ def pprint(di, name_items, fields, str_out=True): items = list(di.keys()) - # ---------- Lenghts - # Lenght of first column + # ---------- Lengths + # Length of first column max_lens = {name_items: max([len(name) for name in items])} - # Lenght of remaining columns + # Length of remaining columns for field in fields: field_vals = [di[name].get(field, '') for name in items] diff --git a/PharmaPy/SimExec.py b/PharmaPy/SimExec.py index ed7826f..0cdf19d 100644 --- a/PharmaPy/SimExec.py +++ b/PharmaPy/SimExec.py @@ -2,7 +2,9 @@ """ Created on Mon Jan 13 12:44:44 2020 -@author: dcasasor +This module contains the SimulationExec class which is responsible for +executing the simulation of a flowsheet in the PharmaPy package. + """ import numpy as np @@ -22,7 +24,61 @@ class SimulationExec: + """ + The SimulationExec class is responsible for executing the simulation of a + flowsheet in the PharmaPy package. + + Parameters + ---------- + pure_path : str + The path to the pure component database. + flowsheet : dict or str + The flowsheet to be simulated. It can be provided as a dictionary or as + a string representation of the flowsheet. + + Attributes + ---------- + NamesSpecies : list + The names of the species in the flowsheet. + StreamTable : pandas.DataFrame + The table containing the stream data. + uos_instances : dict + The instances of the unit operations in the flowsheet. + oper_mode : list + The operating modes of the unit operations. + graph : dict + The graph representation of the flowsheet. + in_degree : dict + The in-degree of each unit operation in the flowsheet. + execution_names : list + The names of the unit operations in the order of execution. + time_processing : dict + The processing time of each unit operation. + result : SimulationResult + The result of the simulation. + connections : dict + The connections between unit operations. + + """ + def __init__(self, pure_path, flowsheet): + """ + Initialize the SimulationExec object. + + Parameters + ---------- + pure_path : str + The path to the pure component database. + flowsheet : dict or str + The flowsheet to be simulated. It can be provided as a dictionary or as + a string representation of the flowsheet. + + Raises + ------ + PharmaPyNonImplementedError + If the provided flowsheet contains recycle stream(s). + + """ # Interfaces thermo_instance = ThermoPhysicalManager(pure_path) @@ -48,6 +104,30 @@ def __init__(self, pure_path, flowsheet): def SolveFlowsheet(self, kwargs_run=None, pick_units=None, verbose=True, steady_state_di=None, tolerances_ss=None, ss_time=0): + """ + Solve the flowsheet simulation. + + Parameters + ---------- + kwargs_run : dict, optional + Additional keyword arguments for running the unit operations. + pick_units : list, optional + The list of unit operations to be solved. If not provided, all unit + operations in the flowsheet will be solved. + verbose : bool, optional + Whether to print verbose output during the simulation. The default is True. + steady_state_di : dict, optional + The dictionary containing steady state information for each unit operation. + tolerances_ss : dict, optional + The dictionary containing tolerances for steady state convergence. + ss_time : float, optional + The steady state time. + + Returns + ------- + None + + """ if kwargs_run is None: kwargs_run = {} @@ -171,7 +251,7 @@ def SetParamEstimation(self, x_data, y_data=None, y_spectra=None, pick_unit=None, **inputs_paramest): """ Set parameter estimation using the aggregated unit operation to a - simulation object + simulation object. Parameters ---------- @@ -208,7 +288,8 @@ def SetParamEstimation(self, x_data, y_data=None, y_spectra=None, control_modifiers : dict, optional Dictionary containing arguments to be passed to a control function. For instance, for a crystallizer with a temperature control with - signature my_control(time, temp_init, ramp): + # instance is already solved, pass data to connection + signature my_control(time, temp_init, ramp): my_modifier = {'temp': {'args': (320, -0.2)}} @@ -313,6 +394,25 @@ def SetParamEstimation(self, x_data, y_data=None, y_spectra=None, def EstimateParams(self, optim_options=None, method='LM', bounds=None, verbose=True): + ''' + Estimate parameters using the parameter estimation object. + + Parameters + ---------- + optim_options : dict, optional + Optimization options. The default is None. + method : str, optional + Optimization method. The default is 'LM'. + bounds : list, optional + Bounds for the optimization. The default is None. + verbose : bool, optional + Whether to print verbose output during the optimization. The default is True. + + Returns + ------- + results : dict + The results of the optimization. + ''' tic = time.time() results = self.ParamInst.optimize_fn(optim_options=optim_options, method=method, @@ -326,8 +426,16 @@ def EstimateParams(self, optim_options=None, method='LM', bounds=None, return results def get_equipment_size(self): - size_equipment = {} + ''' + Function to retrieve the size of the equipment from the simulation data. + Returns + ------- + size_equipment : dict + The size of the equipment in the flowsheet. + ''' + size_equipment = {} + for key, instance in self.uos_instances.items(): if hasattr(instance, 'vol_tot'): size_equipment[key] = instance.vol_tot @@ -342,22 +450,48 @@ def get_equipment_size(self): def GetCAPEX(self, size_equipment=None, k_vals=None, b_vals=None, cepci_vals=None, f_pres=None, f_mat=None, min_capacity=None): + ''' + Function to calculate the CAPEX of the equipment in the flowsheet. + Parameters + ---------- + size_equipment : dict, optional + The size of the equipment in the flowsheet. The default is None. + k_vals : numpy array, optional + The k values for the CAPEX calculation. The default is None. + b_vals : numpy array, optional + The b values for the CAPEX calculation. The default is None. + cepci_vals : numpy array, optional + The CEPCI (Chemical Engineering Plant Cost Index) values for the CAPEX calculation. The default is None. + f_pres : numpy array, optional + The pressure factor for the CAPEX calculation. The default is None. + f_mat : numpy array, optional + The material factor for the CAPEX calculation. The default is None. + min_capacity : numpy array, optional + The minimum capacity for the CAPEX calculation. The default is None. + + Returns + ------- + cost_equip : dict + The CAPEX of the equipment in the flowsheet. + ''' + if size_equipment is None: + # if equipment size is not provided, retrieve from the simulation data size_equipment = self.get_equipment_size() num_equip = len(size_equipment) name_equip = size_equipment.keys() - if cepci_vals is None: + if cepci_vals is None: # cepci stands for Chemical Engineering Plant Cost Index cepci_vals = np.ones(2) - if f_pres is None: + if f_pres is None: # f_pres stands for pressure factor in CAPEX calculations f_pres = np.ones(num_equip) - if f_mat is None: + if f_mat is None: # f_mat stands for material factor in CAPEX calculations f_mat = np.ones(num_equip) - if k_vals is None: + if k_vals is None: # k_vals stands for k1, k2, k3 in CAPEX calculations return size_equipment else: capacities = np.array(list(size_equipment.values())) @@ -366,12 +500,13 @@ def GetCAPEX(self, size_equipment=None, k_vals=None, b_vals=None, a_corr = capacities else: a_corr = np.maximum(min_capacity, capacities) - + # CAPEX calculated by the bare module method k1, k2, k3 = k_vals.T cost_zero = 10**(k1 + k2*np.log10(a_corr) + k3*np.log10(a_corr)**2) b1, b2 = b_vals.T + # bare cost is corrected with material and pressure conditions f_bare = b1 + b2 * f_mat * f_pres cost_equip = cost_zero * f_bare @@ -484,6 +619,8 @@ def get_raw_inlets(self, uo, basis='mass'): dens = stream.getDensity(basis=basis) + total = 0 + if uo.oper_mode == 'Batch': if basis == 'mass': total = stream.mass @@ -549,6 +686,7 @@ def get_raw_inlets(self, uo, basis='mass'): def get_holdup(self, uo, basis='mass'): out = {} + fields = [] if hasattr(uo, '__original_phase__'): phases = uo.__original_phase__ @@ -662,8 +800,13 @@ def GetDuties(self, full_output=False): heat_duties.append(instance.heat_duty) equipment_ids.append(key) - heat_duties = np.array(heat_duties) - heat_duties = pd.DataFrame(heat_duties, index=equipment_ids, + # heat_duties = np.array(heat_duties) + if heat_duties: + heat_duties = np.array(heat_duties) + heat_duties_df = pd.DataFrame(heat_duties, index=equipment_ids, columns=['heating', 'cooling']) + else: + heat_duties_df = pd.DataFrame(columns=['heating', 'cooling']) + heat_duties = pd.DataFrame(heat_duties, index=equipment_ids, columns=['heating', 'cooling']) duties_ids = np.array(duty_ids) @@ -697,14 +840,25 @@ def GetOPEX(self, cost_raw, include_holdups=True, steady_raw=False, duty_unit_cost = np.zeros_like(map_duties, dtype=np.float64) for ind, row in enumerate(map_duties): duty_unit_cost[ind] = heat_exchange_cost[row] - - duty_cost = np.abs(duties)*1e-9 * duty_unit_cost + + # duty_cost = np.abs(duties)*1e-9 * duty_unit_cost + if duties.size > 0: + duty_cost = np.abs(duties) * 1e-9 * duty_unit_cost + else: + duty_cost = np.array([]) # ---------- Raw materials - raw_materials = self.GetRawMaterials( - include_holdups, steady_raw, **kwargs_items.get('raw_materials', - {})) - raw_cost = cost_raw * raw_materials + raw_materials = self.GetRawMaterials(basis = 'mass', totals=include_holdups, steady_state = steady_raw) + # raw_materials = self.GetRawMaterials( + # include_holdups, steady_raw, **kwargs_items.get('raw_materials', + # {})) + + # raw_materials is already cutting off non-mass fields in get_holdup/get_raw_inlets + selected_cols = raw_materials.columns[1:] + if raw_materials.size > 0: + raw_cost = cost_raw * raw_materials[selected_cols] + else: + raw_cost = pd.DataFrame(columns=selected_cols) # ---------- Labor labor_cost = self.GetLabor(**kwargs_items.get('labor', {})) diff --git a/install_instructions.txt b/install_instructions.txt index be19104..cfc2ceb 100644 --- a/install_instructions.txt +++ b/install_instructions.txt @@ -2,7 +2,7 @@ After installing miniconda, do the following: Step 1: - WINDOWS: Open an Anaconda PowerShell ("powershell" for short) by clicking on: Start --> Anaconda3 (XX-bit) --> Anaconda PowerShell Prompt (Miniconda3). Right clic on it and execute as an administrator + WINDOWS: Open an Anaconda PowerShell ("powershell" for short) by clicking on: Start --> Anaconda3 (XX-bit) --> Anaconda PowerShell Prompt (Miniconda3). Right click on it and execute as an administrator LINUX/MAC: Open a new terminal where the this file is. Step 2 (WINDOWS only): Copy (CTRL+C) the path to the installation files (where this file is). This can be done by clicking on the top horizontal bar of the file explorer, which will highlight the path we need to go to. For instance, after copying, my path is C:\Users\dcasasor\OneDrive - purdue.edu\postdoc\PharmaPy_dcasasor @@ -10,10 +10,11 @@ Step 2 (WINDOWS only): Copy (CTRL+C) the path to the installation files (where t Step 3 (WINDOWS only): Introduce the following command in the powershell: cd '' and hit Enter. Do not include the <> characters, but do include the quotes. Replace text in between the <> characters with the path you copied in Step 2. Use CTRL+V to avoid transcription errors Step 4: After this step, both PharmaPy dependencies and PharmaPy itself should be installed - WINDOWS: Do .\InstallOnWindows.bat and hit Enter. Follow the instructions on screen - LINUX/MAC: To give permissions to the installation file, do chmod +x InstallOnMac.sh. Then, execute it with ./InstallOnMac.sh. Follow the instructions on screen + WINDOWS: Do .\InstallOnWindows.bat and hit Enter. Follow the instructions on screen + LINUX/MAC: To give permissions to the installation file, do chmod +x InstallOnMac.sh. + Then, execute it with source InstallOnMac.sh. Follow the instructions on screen -Step 5: Activate the new environment bydoing: conda activate (exclude <>) +Step 5: Activate the new environment by doing: conda activate (exclude <>) Step 6: Test PharmaPy by moving to the tests/ directory (cd tests) and then running: python reactor_tests.py @@ -39,7 +40,7 @@ the analysis we need to do! If you want to use base, you may continue to step When user does not have a write permission to a required path, please direct the path of the directory -step 1.1.1: Disclose the available enviornment directories. +step 1.1.1: Disclose the available environment directories. conda config -- show envs_dirs diff --git a/requirements.txt b/requirements.txt index ce3a782..f23b2f6 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,7 +1,7 @@ numpy +cython < 3 scipy matplotlib pandas autograd assimulo -cyipopt diff --git a/setup.py b/setup.py index bdc2fa0..9eceef7 100644 --- a/setup.py +++ b/setup.py @@ -1,22 +1,15 @@ -from setuptools import setup, find_packages +from setuptools import setup, find_packages, Extension -setuptools_kwargs = { - 'zip_safe': False, - 'install_requires': ['pyomo', - 'coverage', - 'numpy', - 'scipy', - 'pandas', - 'assimulo'], - 'scripts': [], - 'include_package_data': True -} +# Read the content of requirements.txt into a list +with open("requirements.txt", "r") as f: + requirements = f.read().splitlines() setup(name='PharmaPy', - version='0.0.0', + version='0.0.1', packages=find_packages(), author='Daniel Casas-Orozco', author_email='dcasasor@purdue.edu', license='', url='', - **setuptools_kwargs) + py_modules=["PharmaPy"], + install_requires=requirements) \ No newline at end of file