diff --git a/biorefineries/microalgae/TRY_analysis.py b/biorefineries/microalgae/TRY_analysis.py new file mode 100644 index 00000000..368f2312 --- /dev/null +++ b/biorefineries/microalgae/TRY_analysis.py @@ -0,0 +1,162 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Titer, Rate, Yield (TRY) analysis for microalgae biorefinery + +Based on succinic project but adapted for microalgae system structure +""" + +import numpy as np +import pandas as pd +import os +from .uncertainties import create_model + +def run_TRY_analysis(): + """Run Titer, Rate, Yield analysis for microalgae system""" + print("Creating model for TRY analysis...") + model, namespace_dict = create_model() + + # Load parameter distributions + current_dir = os.path.dirname(__file__) + param_file = os.path.join(current_dir, 'parameter_distributions.xlsx') + + if not os.path.exists(param_file): + raise FileNotFoundError(f"Parameter distributions file not found: {param_file}") + + print(f"Loading parameter distributions from {param_file}...") + model.load_parameter_distributions(param_file, namespace_dict) + + # Define TRY parameter ranges based on microalgae system + # These ranges should be adjusted based on your specific system + titer_range = np.linspace(1.0, 4.0, 11) # g/L - caproic acid titer + yield_range = np.linspace(0.15, 0.35, 11) # g/g - C6 yield + + results = [] + + # Find TRY parameters from loaded parameters + titer_param = None + yield_param = None + + for param in model.parameters: + param_name_lower = param.name.lower() + if 'titer' in param_name_lower and 'c6' in param_name_lower: + titer_param = param + print(f"Found titer parameter: {param.name}") + elif 'yield' in param_name_lower and 'c6' in param_name_lower: + yield_param = param + print(f"Found yield parameter: {param.name}") + + if titer_param is None: + print("Warning: Could not find titer parameter. Available parameters:") + for param in model.parameters: + print(f" - {param.name}") + # Use first parameter with 'titer' in name as fallback + for param in model.parameters: + if 'titer' in param.name.lower(): + titer_param = param + print(f"Using fallback titer parameter: {param.name}") + break + + if yield_param is None: + print("Warning: Could not find yield parameter. Available parameters:") + for param in model.parameters: + print(f" - {param.name}") + # Use first parameter with 'yield' in name as fallback + for param in model.parameters: + if 'yield' in param.name.lower(): + yield_param = param + print(f"Using fallback yield parameter: {param.name}") + break + + if titer_param is None or yield_param is None: + print("Error: Could not find both titer and yield parameters") + return None + + # Run TRY analysis + print(f"\nRunning TRY analysis...") + print(f"Titer range: {min(titer_range):.2f} - {max(titer_range):.2f}") + print(f"Yield range: {min(yield_range):.3f} - {max(yield_range):.3f}") + + baseline_sample = model.get_baseline_sample() + total_combinations = len(titer_range) * len(yield_range) + completed = 0 + + for i, titer in enumerate(titer_range): + for j, yield_val in enumerate(yield_range): + sample = baseline_sample.copy() + sample[titer_param.index] = titer + sample[yield_param.index] = yield_val + + try: + result = model(sample) + result['Titer'] = titer + result['Yield'] = yield_val + results.append(result) + completed += 1 + + if completed % 10 == 0: + print(f"Completed {completed}/{total_combinations} combinations") + + except Exception as e: + print(f"Error at titer={titer:.2f}, yield={yield_val:.3f}: {e}") + continue + + if results: + results_df = pd.DataFrame(results) + print(f"\nTRY analysis completed! {len(results)} successful evaluations out of {total_combinations} total combinations") + + # Print summary statistics + print("\nSummary statistics:") + for col in ['MFSP', 'GWP', 'FEC', 'TCI']: + if col in results_df.columns: + print(f"{col}: min={results_df[col].min():.4f}, max={results_df[col].max():.4f}, mean={results_df[col].mean():.4f}") + + return results_df + else: + print("No successful evaluations in TRY analysis") + return None + +def analyze_TRY_results(results_df, save_path=None): + """Analyze TRY results and identify optimal regions""" + if results_df is None: + print("No results to analyze") + return + + print("\nAnalyzing TRY results...") + + # Find optimal conditions (minimum MFSP) + if 'MFSP' in results_df.columns: + min_mfsp_idx = results_df['MFSP'].idxmin() + optimal_conditions = results_df.loc[min_mfsp_idx] + print(f"\nOptimal conditions (minimum MFSP):") + print(f" Titer: {optimal_conditions['Titer']:.2f}") + print(f" Yield: {optimal_conditions['Yield']:.3f}") + print(f" MFSP: {optimal_conditions['MFSP']:.4f} $/kg") + if 'GWP' in results_df.columns: + print(f" GWP: {optimal_conditions['GWP']:.4f} kg CO2-eq/kg") + + # Save results if path provided + if save_path: + results_df.to_csv(save_path, index=False) + print(f"\nResults saved to: {save_path}") + +if __name__ == '__main__': + try: + results = run_TRY_analysis() + + if results is not None: + # Analyze results + analyze_TRY_results(results) + + # Save results + results_dir = 'results' + if not os.path.exists(results_dir): + os.makedirs(results_dir) + + save_path = os.path.join(results_dir, 'TRY_analysis_results.csv') + analyze_TRY_results(results, save_path) + + except Exception as e: + print(f"Error in TRY analysis: {e}") + import traceback + traceback.print_exc() \ No newline at end of file diff --git a/biorefineries/microalgae/_chemicals.py b/biorefineries/microalgae/_chemicals.py new file mode 100644 index 00000000..3616db2d --- /dev/null +++ b/biorefineries/microalgae/_chemicals.py @@ -0,0 +1,405 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Sat June 21 20:50:00 2025 + +Microalgae biorefinery to produce medium chain fatty acids +by anaerobic fermentation without external electron donor addition- chemicals database + +References +---------- +[1] BioSTEAM Documentation: + https://biosteam.readthedocs.io/en/latest/API/thermosteam/Chemicals.html +[2] Cortes-Peña et al., BioSTEAM: A Fast and Flexible Platform for the Design, + Simulation, and Techno-Economic Analysis of Biorefineries under Uncertainty. + ACS Sustainable Chem. Eng. 2020, 8 (8), 3302–3310. +[3] 3-Hydroxypropionic acid biorefineries project: + https://github.com/BioSTEAMDevelopmentGroup/Bioindustrial-Park/tree/master/biorefineries/HP + +@author: Xingdong Shi +@version: 0.0.2 +""" + +#import biosteam as bst +import thermosteam as tmo +from thermosteam import functional as fn +#from biorefineries.cellulosic import chemicals as cellulosic_chems +#import biorefineries.sugarcane as sc +#from biorefineries.sugarcane import chemicals as sugarcane_chems +from biorefineries.cornstover import create_chemicals as cs_create_chemicals +from biorefineries.corn import create_chemicals as corn_create_chemicals + +# Import chemicals from other biorefineries +cornstover_chems = cs_create_chemicals() +corn_chems = corn_create_chemicals() + +# Create microalgae chemicals system +chems = Microalgae_chemicals = tmo.Chemicals([]) + +# Track chemicals by creation method +database_chemicals_dict = {} +copied_chemicals_dict = {} +defined_chemicals_dict = {} + +def chemical_database(ID, phase=None, **kwargs): + """Create chemical from database""" + chemical = tmo.Chemical(ID, **kwargs) + if phase: + chemical.at_state(phase) + chemical.phase_ref = phase + chems.append(chemical) + database_chemicals_dict[ID] = f'{ID}: {chemical.formula}/{chemical.MW}' + return chemical + + +def chemical_copied(ID, ref_chemical, **data): + """Copy existing chemical""" + chemical = ref_chemical.copy(ID) + chems.append(chemical) + for i, j in data.items(): + setattr(chemical, i, j) + copied_chemicals_dict[ID] = f'{ID}: {chemical.formula}/{chemical.MW}' + return chemical + +def chemical_defined(ID, **kwargs): + """Define custom chemical""" + chemical = tmo.Chemical.blank(ID, **kwargs) + chems.append(chemical) + defined_chemicals_dict[ID] = f'{ID}: {chemical.formula}/{chemical.MW}' + return chemical + +# ============================================================================= +# Constants +# ============================================================================= + +# Microalgae properties based on Nannochloropsis salina +# https://www.sciencedirect.com/science/article/pii/S0378382016302983 +Cp = 3.949 # J/g/K +rho = 1.0093 # g/cm3 +_cal2joule = 4.184 + +# ============================================================================= +# Basic chemicals +# ============================================================================= + +H2O = chemical_database('H2O') + +# ============================================================================= +# Gases +# ============================================================================= + +O2 = chemical_database('O2', phase='g', Hf=0) +N2 = chemical_database('N2', phase='g', Hf=0) +CH4 = chemical_database('CH4', phase='g') +CO2 = chemical_database('CO2', phase='g') +H2 = chemical_database('H2', phase='g', Hf=0) +NH3 = chemical_database('NH3', phase='g', Hf=-10963*_cal2joule) +NitricOxide = chemical_database('NitricOxide', phase='g') +NO2 = chemical_database('NO2', phase='g') +H2S = chemical_database('H2S', phase='g', Hf=-4927*_cal2joule) +SO2 = chemical_database('SO2', phase='g') + +# ============================================================================= +# Solids for combustion +# ============================================================================= + +P4O10 = chemical_database('P4O10', phase='s', Hf=-713.2*_cal2joule) + +# ============================================================================= +# Soluble inorganics +# ============================================================================= + +HCl = chemical_database('HCl', phase='l') +H2SO4 = chemical_database('H2SO4', phase='l') +HNO3 = chemical_database('HNO3', phase='l', Hf=-41406*_cal2joule) +NaOH = chemical_database('NaOH', phase='l') +NH4OH = chemical_database('NH4OH', phase='l') + +# Lime (Calcium hydroxide) for pH control and boiler +Lime = chemical_database('Lime', search_ID='CalciumDihydroxide', phase='l') + +# Calcium sulfate for boiler ash +CaSO4 = chemical_database('CaSO4', phase='s', Hf=-342531*_cal2joule) +# Use Lastovka solid model instead of default Perry 151 model +CaSO4.Cn.move_up_model_priority('LASTOVKA_S', 0) + +# Sodium sulfate for neutralization reactions +Na2SO4 = chemical_database('Na2SO4', phase='l', Hf=-1356.38e3) +AmmoniumSulfate = chemical_database('AmmoniumSulfate', phase='l', Hf=-1186.8e3) # 新增硫酸铵，主流BioSTEAM项目标准 + +# ============================================================================= +# Products +# ============================================================================= + +# Fatty acids +AceticAcid = chemical_database('AceticAcid') +PropionicAcid = chemical_database('PropionicAcid') +ButyricAcid = chemical_database('ButyricAcid') +ValericAcid = chemical_database('ValericAcid') +CaproicAcid = chemical_database('CaproicAcid') +HeptanoicAcid = chemical_database('HeptanoicAcid') +CaprylicAcid = chemical_database('CaprylicAcid') + +# Alcohols +Ethanol = chemical_database('Ethanol') +Butanol = chemical_database('Butanol') +Hexanol = chemical_database('Hexanol') + + +# Sugars and carbohydrates +Glucose = chemical_database('Glucose', phase='l') +Fructose = chemical_database('Fructose', phase='l') + +# Sugar oligomers +GlucoseOligomer = chemical_defined('GlucoseOligomer', phase='l', formula='C6H10O5', + Hf=-233200*_cal2joule) +GlucoseOligomer.copy_models_from(Glucose, ['Hvap', 'Psat', 'Cn', 'mu', 'kappa']) +Extract = chemical_copied('Extract', Glucose) + +Xylose = chemical_database('Xylose') +Xylose.copy_models_from(Glucose, ['Hvap', 'Psat', 'mu']) +XyloseOligomer = chemical_defined('XyloseOligomer', phase='l', formula='C5H8O4', + Hf=-182100*_cal2joule) +XyloseOligomer.copy_models_from(Xylose, ['Hvap', 'Psat', 'Cn', 'mu']) + +Sucrose = chemical_database('Sucrose', phase='l') +Sucrose.Cn.move_up_model_priority('DADGOSTAR_SHAW', 0) +Cellobiose = chemical_database('Cellobiose', phase='l', Hf=-480900*_cal2joule) + +Mannose = chemical_database('Mannose', phase='l', Hf=Glucose.Hf) +Mannose.copy_models_from(Glucose, ['Hvap', 'Psat', 'Cn', 'mu']) +MannoseOligomer = chemical_copied('MannoseOligomer', GlucoseOligomer) + +Galactose = chemical_database('Galactose', phase='l', Hf=Glucose.Hf) +Galactose.copy_models_from(Glucose, ['Hvap', 'Psat', 'Cn','mu']) +GalactoseOligomer = chemical_copied('GalactoseOligomer', GlucoseOligomer) + +Arabinose = chemical_database('Arabinose', phase='l', Hf=Xylose.Hf) +Arabinose.copy_models_from(Xylose, ['Hvap', 'Psat', 'mu']) +ArabinoseOligomer = chemical_copied('ArabinoseOligomer', XyloseOligomer) + +SolubleLignin = chemical_database('SolubleLignin', search_ID='Vanillin', + phase='l', Hf=-108248*_cal2joule) + +# ============================================================================= +# Extractants and solvents +# ============================================================================= +Hexane = chemical_database('Hexane', phase='l') +Heptane = chemical_database('Heptane', phase='l') +Octane = chemical_database('Octane', phase='l') +Toluene = chemical_database('Toluene', phase='l') +Benzene = chemical_database('Benzene', phase='l') +Chloroform = chemical_database('Chloroform', phase='l') +Dichloromethane = chemical_database('Dichloromethane', phase='l') +Octanol = chemical_database('Octanol', phase= 'l') +OleylAlcohol = chemical_database('OleylAlcohol', phase= 'l') + +# Ionic liquid extractants +TOA = chemical_database('TOA', search_ID='tri-n-octylamine') +AQ336 = chemical_database('AQ336', search_ID='63393-96-4') # aliquat 336 +AQ336.copy_models_from(TOA, ('Psat', 'Hvap', 'V')) +AQ336._Dortmund = TOA.Dortmund +AQ336.Hfus = TOA.Hfus + +# ============================================================================= +# Microalgae components +# ============================================================================= + +# Protein - main microalgae component +Protein = chemical_defined('Protein', + phase='s', + formula='CH1.57O0.31N0.29S0.007', + Hf=-17618*_cal2joule) + +# Carbohydrate +Carbohydrate = chemical_defined('Carbohydrate', + phase='s', + formula='CH2O', + MW=180.16, + Hf=-233000) + +# Ash +Ash = chemical_defined('Ash', + phase='s', + MW=1, + Hf=0, + HHV=0, + LHV=0) + + +WWTsludge = chemical_defined('WWTsludge', phase='s', formula='CH1.66O0.33N0.22', MW=24.6, Hf=-110000) +# ============================================================================= +# Enzymes and microorganisms +# ============================================================================= + +# Base enzyme +BaseEnzyme = chemical_defined('BaseEnzyme', + phase='s', + formula='CH1.57O0.31N0.29S0.007', + Hf=-17618*_cal2joule, + HHV=0, + LHV=0) + +# Specific enzymes +AlphaAmylase = chemical_copied('AlphaAmylase', BaseEnzyme, Hf=-2000000) +GlucoAmylase = chemical_copied('GlucoAmylase', BaseEnzyme, Hf=-2500000) +Lipase = chemical_copied('Lipase', BaseEnzyme, Hf=-1800000) +Protease = chemical_copied('Protease', BaseEnzyme, Hf=-2200000) + +# Fermentation microorganisms +FermMicrobe = chemical_defined('FermMicrobe', phase='l', + formula='CH1.78O0.44N0.24', Hf=-103310.) + +# ============================================================================= +# Insoluble organics +# ============================================================================= + +Glucan = chemical_defined('Glucan', phase='s', formula='C6H10O5', Hf=-233200*_cal2joule) +Glucan.copy_models_from(Glucose, ['Cn']) +Mannan = chemical_copied('Mannan', Glucan) +Galactan = chemical_copied('Galactan', Glucan) + +Xylan = chemical_defined('Xylan', phase='s', formula='C5H8O4', Hf=-182100*_cal2joule) +Xylan.copy_models_from(Xylose, ['Cn']) +Arabinan = chemical_copied('Arabinan', Xylan) + +Cellulose = chemical_database('Cellulose', phase='s') +Hemicellulose = chemical_database('Hemicellulose', phase='s') + +# ============================================================================= +# Mixtures +# ============================================================================= +# Boiler chemicals +BoilerChems = chemical_database('BoilerChems', search_ID='DiammoniumPhosphate', + phase='l', Hf=0, HHV=0, LHV=0) + +# ============================================================================= +# Fillers +# ============================================================================= + +BaghouseBag = chemical_defined('BaghouseBag', phase='s', MW=1, Hf=0, HHV=0, LHV=0) +BaghouseBag.Cn.add_model(0) +CoolingTowerChems = chemical_copied('CoolingTowerChems', BaghouseBag) + +Microalgae = chemical_defined( + 'Microalgae', + phase='s', + formula='C3.2H0.4O1.5N0.3', # comsuption from microalgal biomass componets https://www.sciencedirect.com/science/article/pii/S0196890419313184?via%3Dihub + Hf=-2138000, +) + + + + +# ============================================================================= +# Chemical groups for simulation +# ============================================================================= + +chemical_groups = dict( + OtherSugars = ('Arabinose', 'Mannose', 'Galactose', 'Cellobiose', + 'Sucrose', 'Fructose'), + SugarOligomers = ('GlucoseOligomer', 'XyloseOligomer', 'GalactoseOligomer', + 'ArabinoseOligomer', 'MannoseOligomer'), + InorganicSolubleSolids = ('NaOH', 'HNO3', 'NaNO3', 'BoilerChems', + 'Na2SO4', 'AmmoniumSulfate'), + Furfurals = ('Furfural', 'HMF'), + OtherOrganics = ('Xylitol', 'Glycerol'), + COSOxNOxH2S = ('NitricOxide', 'NO2', 'SO2', 'H2S'), + Proteins = ('Protein', 'BaseEnzyme', 'AlphaAmylase', + 'GlucoAmylase', 'Lipase', 'Protease', 'Microalgae'), + CellMass = ('FermMicrobe', 'WWTsludge'), + OtherInsolubleSolids = ('Ash', 'CalciumDihydroxide', 'CaSO4', + 'BaghouseBag', 'CoolingTowerChems'), + OtherStructuralCarbohydrates = ('Glucan', 'Xylan', 'Arabinan', + 'Mannan', 'Galactan', 'Cellulose', 'Hemicellulose'), + SeparatelyListedOrganics = ('Ethanol', 'Glucose', 'Xylose', 'AceticAcid'), + SpearatedlyListedOthers = ('H2O', 'NH3', 'H2SO4', 'CO2', 'CH4', 'O2', 'N2'), + MCCA = ('CaproicAcid', 'HeptanoicAcid', 'CaprylicAcid'), + SCFA = ('AceticAcid', 'PropionicAcid', 'ButyricAcid', 'ValericAcid'), + Alcohols = ('Ethanol', 'Butanol', 'Hexanol'), + Extractants = ('Hexane', 'Heptane', 'Octane', 'Toluene', + 'Benzene', 'Chloroform', 'Dichloromethane', + 'TOA', 'AQ336', 'Octanal') +) +phase_change_chemicals = ['Methanol', 'MCCA', 'SCFA', 'H2O', 'MethylAcetate', 'Denaturant', + 'AceticAcid', 'MethylAcetate', 'MethylLactate', + 'EthylLactate', 'Furfural', 'MethylSuccinate', + 'SuccinicAcid', 'LacticAcid', 'HMF', 'Alcohols', 'Extractants'] + + +soluble_groups = ('OtherSugars', 'SugarOligomers', + 'Furfurals', 'OtherOrganics', 'Proteins', 'CellMass', + 'SeparatelyListedOrganics') +soluble_organics = list(sum([chemical_groups[i] for i in soluble_groups], ())) + +solubles = tuple(soluble_organics) + chemical_groups['InorganicSolubleSolids'] + ('H2SO4',) + +insoluble_groups = ('OtherInsolubleSolids', 'OtherStructuralCarbohydrates') +insolubles = sum([chemical_groups[i] for i in insoluble_groups], ('WWTsludge',)) + +for chem in chems: + if chem.ID in phase_change_chemicals: pass + elif chem.locked_state: pass + else: + # Set phase_ref to avoid missing model errors + if chem.phase_ref == 'g': + chem.at_state('g') + if chem.ID in solubles: + chem.phase_ref = 'l' + chem.at_state('l') + if chem.ID in insolubles: + chem.phase_ref = 's' + chem.at_state('s') + +# Set chemical heat capacity +# Cp of biomass (1.25 J/g/K) from Leow et al., Green Chemistry 2015, 17 (6), 3584–3599 +for chemical in (Protein, FermMicrobe, BaseEnzyme): + chemical.Cn.add_model(1.25*chemical.MW) + +# Set chemical molar volume following assumptions in lipidcane biorefinery, +# assume densities for solulables and insolubles to be 1e5 and 1540 kg/m3, respectively + +def set_rho(chemical, rho): + V = fn.rho_to_V(rho, chemical.MW) + chemical.V.add_model(V, top_priority=True) + +for chemical in chems: + if chemical.ID in phase_change_chemicals: pass + elif chemical.ID in solubles: set_rho(chemical, 1e5) + elif chemical.ID in insolubles: set_rho(chemical, 1540) + +# The Lakshmi Prasad model gives negative kappa values for some chemicals +for chemical in chems: + if chemical.locked_state: + try: chemical.kappa.move_up_model_priority('Lakshmi Prasad', -1) + except: pass + +# Default missing properties of chemicals to those of water, +for chemical in chems: chemical.default() + + +# ============================================================================= +# Set Hfus to 0 if None +# ============================================================================= + +for chem in Microalgae_chemicals: + if chem.Hfus is None: + chem.Hfus = 0 + +# ============================================================================= +# Add chemicals from other biorefineries +# ============================================================================= + +chems.append(corn_chems.Starch) +chems.append(corn_chems.Fiber) +chems.append(corn_chems.SolubleProtein) +chems.append(corn_chems.InsolubleProtein) +chems.append(corn_chems.Oil) +chems.append(corn_chems.Yeast) + + +chems.compile() +tmo.settings.set_thermo(chems) + + + diff --git a/biorefineries/microalgae/analysis_utils.py b/biorefineries/microalgae/analysis_utils.py new file mode 100644 index 00000000..6c311e42 --- /dev/null +++ b/biorefineries/microalgae/analysis_utils.py @@ -0,0 +1,110 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Microalgae biorefinery analysis utilities + +Helper functions for creating LCA objects and other analysis utilities. + +@author: Xingdong Shi +@version: 0.0.1 +""" + +from . import lca + +def create_microalgae_lca_simple(system, tea=None): + """ + Simplified function to create microalgae LCA object with default parameters. + + Parameters + ---------- + system : System + Microalgae production system + tea : TEA, optional + TEA object (not used in current implementation) + + Returns + ------- + LCA + Life cycle assessment object + """ + + # Find main product + main_product = None + for stream in system.products: + if stream.get_total_flow('kg/hr') > 0: + main_product = stream + break + + if main_product is None: + # If no products found, try to find caproic acid stream + for stream in system.streams: + if 'caproic' in stream.ID.lower(): + main_product = stream + break + + # Find boiler + boiler = None + for unit in system.units: + if 'BT' in unit.ID or 'boiler' in unit.ID.lower(): + boiler = unit + break + + # Default main product chemical IDs + main_product_chemical_IDs = ['CaproicAcid'] + + # Create LCA object + microalgae_lca = lca.create_microalgae_lca( + system=system, + main_product=main_product, + main_product_chemical_IDs=main_product_chemical_IDs, + boiler=boiler + ) + + return microalgae_lca + +def get_fermentation_unit(system): + """ + Find the fermentation unit in the system. + + Parameters + ---------- + system : System + Microalgae production system + + Returns + ------- + Unit or None + Fermentation unit if found, None otherwise + """ + for unit in system.units: + if ('ferment' in unit.ID.lower() or + 'MCCA' in unit.ID or + 'R' in unit.ID and hasattr(unit, 'conversion')): + return unit + return None + +def get_main_product_stream(system): + """ + Find the main product stream in the system. + + Parameters + ---------- + system : System + Microalgae production system + + Returns + ------- + Stream or None + Main product stream if found, None otherwise + """ + # First try to find streams with significant flow + for stream in system.products: + if stream.get_total_flow('kg/hr') > 0: + return stream + + # If no products found, try to find caproic acid stream + for stream in system.streams: + if 'caproic' in stream.ID.lower() and stream.get_total_flow('kg/hr') > 0: + return stream + + return None \ No newline at end of file diff --git a/biorefineries/microalgae/facilities.py b/biorefineries/microalgae/facilities.py new file mode 100644 index 00000000..92e1eb03 --- /dev/null +++ b/biorefineries/microalgae/facilities.py @@ -0,0 +1,646 @@ +# -*- coding: utf-8 -*- +""" +Created on Sat June 21 10:50:00 2025 + +Microalgae biorefinery to produce medium chain fatty acids +by anaerobic fermentation without external electron donor addition- facilities + +References +---------- +[1] BioSTEAM Documentation: + https://biosteam.readthedocs.io/en/latest/tutorial/Creating_a_System.html +[2] Cortes-Peña et al., BioSTEAM: A Fast and Flexible Platform for the Design, + Simulation, and Techno-Economic Analysis of Biorefineries under Uncertainty. + ACS Sustainable Chem. Eng. 2020, 8 (8), 3302–3310. +[3] 3-Hydroxypropionic acid biorefineries project: + https://github.com/BioSTEAMDevelopmentGroup/Bioindustrial-Park/tree/master/biorefineries/HP +[4] Succinic projest + https://github.com/BioSTEAMDevelopmentGroup/Bioindustrial-Park/tree/master/biorefineries/succinic + +@author: Xingdong Shi +@version: 0.0.1 +""" + +import biosteam as bst +import thermosteam as tmo +from biosteam import HeatUtility, Facility +from biosteam.units.decorators import cost +from thermosteam import Stream +from .utils import CEPCI +import flexsolve as flx + +__all__ = ('CIP', 'ADP', 'CT', 'PWC', 'BT', 'HXNWithMin') + + +# %% + +# ============================================================================= +# Clean-in-place system +# ============================================================================= +@cost(basis='Flow rate', ID='System', units='kg/hr', + cost=421000, S=63, CE=CEPCI[2009], n=0.6, BM=1.8) +class CIP(Facility): + network_priority = 3 + line = 'Clean-in-place system' + + +# %% + +# ============================================================================= +# Air distribution package, size based on stream 950 in Humbird et al. +# ============================================================================= +@cost(basis='Flow rate', ID='Plant air compressor', units='kg/hr', + kW=111.855, cost=28000, S=83333, CE=CEPCI[2010], n=0.6, BM=1.6) +@cost(basis='Flow rate', ID='Plant air reciever', units='kg/hr', + cost=16000, S=83333, CE=CEPCI[2009], n=0.6, BM=3.1) +@cost(basis='Flow rate', ID='Instrument air dryer', units='kg/hr', + cost=15000, S=83333, CE=CEPCI[2009], n=0.6, BM=1.8) +class ADP(Facility): + network_priority = 3 + line = 'Air distribution package' + + def __init__(self, ID='', ins=None, outs=(), ratio=None): + Facility.__init__(self, ID, ins, outs) + self.ratio = ratio + + def _design(self): + self.design_results['Flow rate'] = 83333 * self.ratio + + +# %% + +# ============================================================================= +# Chilled water package +# ============================================================================= + +# Use BioSTEAM built-in ChilledWaterPackage class +CWP = bst.ChilledWaterPackage + + +# %% + +# ============================================================================= +# Cooling tower +# ============================================================================= + +# Use BioSTEAM built-in CoolingTower class +CT = bst.CoolingTower + + +# %% + +# ============================================================================= +# Process water center +# ============================================================================= + +@cost(basis='Flow rate', ID='Tank', units='kg/hr', + cost=250000, S=451555, CE=CEPCI[2009], n=0.7, BM=1.7) +@cost(basis='Flow rate', ID='Circulating pump', units='kg/hr', + kW=55.9275, cost=15292, S=518924, CE=CEPCI[2010], n=0.8, BM=3.1) +@cost(basis='Flow rate', ID='Makeup water pump', units='kg/hr', + kW=14.914, cost=6864, S=155564, CE=CEPCI[2010], n=0.8, BM=3.1) +class PWC(Facility): + _N_ins = 2 + _N_outs = 2 + _units= {'Flow rate': 'kg/hr'} + + network_priority = 2 + line = 'Process water center' + + def __init__(self, ID='', ins=None, outs=(), process_water_streams=None, + recycled_blowdown_streams=None): + Facility.__init__(self, ID, ins, outs) + self.process_water_streams = process_water_streams + self.recycled_blowdown_streams = recycled_blowdown_streams + + def _run(self): + makeup, RO_water = self.ins + process_water, discharged = self.outs + + water_demand = sum(i.imol['Water'] for i in self.process_water_streams) + water_needs = water_demand - RO_water.imol['Water'] + self.recycled_water = RO_water.imass['Water'] + + if self.recycled_blowdown_streams: + water_needs -= sum(i.imol['Water'] for i in self.recycled_blowdown_streams) + self.recycled_water += sum(i.imass['Water'] for i in self.recycled_blowdown_streams) + + if water_needs > 0: + makeup.imol['Water'] = water_needs + discharged.empty() + else: + discharged.imol['Water'] = - water_needs + makeup.empty() + + process_water.mol = makeup.mol + RO_water.mol - discharged.mol + + self.design_results['Flow rate'] = self.F_mass_in + +# ============================================================================= +# Boiler and turbogenerator +# ============================================================================= + + +@cost(basis='Flow rate', ID='Boiler', units='kg/hr', + cost=36500000, S=239000, CE=CEPCI[2009], n=0.8, BM=1.8) +class BT(bst.Facility): + """ + Create a BoilerTurbogenerator object that will calculate electricity + generation from burning the feed. It also takes into account how much + steam is being produced, and the required cooling utility of the turbo + generator. All capital cost correlations are based on [1]_. + + Parameters + ---------- + ins : + * [0] Liquid/solid feed that will be burned. + * [1] Gas feed that will be burned. + * [2] Make-up water. + * [3] Natural gas/fuel to satisfy steam and electricity demand. + * [4] Lime for flue gas desulfurization. + * [5] Boiler chemicals. + * [6] Air or oxygen-rich gas. + outs : + * [0] Total emissions produced. + * [1] Blowdown water. + * [2] Ash disposal. + boiler_efficiency : float, optional + Fraction of heat transferred to steam. Defaults to 0.8. + turbo_generator_efficiency : float, optional + Fraction of steam heat converted to electricity. Defaults to 0.85. + agent : UtilityAgent, optional + Steam produced. Defaults to low pressure steam. + other_agents = () : Iterable[UtilityAgent], optional + Other steams produced. Defaults to all other heating agents. + fuel_source : str, optional + Name fuel used to satisfy steam and electricity demand. Defaults to 'CH4'. + fuel_price : float, optional + Price of natural gas [USD/kg]. Same as `bst.stream_utility_prices['Natural gas']`, + which defaults to 0.218. + ash_disposal_price : float, optional + Price of disposing ash [USD/kg]. Same as `bst.stream_utility_prices['Ash disposal']`, + which defaults to -0.0318. + satisfy_system_electricity_demand : bool, optional + Whether to purchase natural gas to satisfy system electricity demand + if there is not enough heat from process feeds (i.e., inlets 0 and 1). + If True, natural gas is purchased to satisfy system heat and electricity demand + when there is not enough heat from the feed and gas. + If False, natural gas is only purchased to satisfy system heat demand + and electricity will be purchased from the grid if there is not + enough heat from the feeds. + In either case, if there is excess heat from the process feeds, + electricity will still be produced. + boiler_efficiency_basis : str, optional + Basis of boiler efficiency. Defaults to 'LHV' (i.e., lower heating value). + 'HHV' (i.e., higher heating value) is also a valid basis. + + Examples + -------- + Create a boiler-turbogenerator system that uses sugarcane bagasse to + produce steam for a distillation unit and any excess steam for surplus electricity: + + >>> import biosteam as bst + >>> from biorefineries import cane + >>> chemicals = cane.create_sugarcane_chemicals() + >>> chemicals.define_group( + ... name='Fiber', + ... IDs=['Cellulose', 'Hemicellulose', 'Lignin'], + ... composition=[0.4704 , 0.2775, 0.2520], + ... wt=True, # Composition is given as weight + ... ) + >>> bst.settings.set_thermo(chemicals) + >>> dilute_ethanol = bst.Stream('dilute_ethanol', Water=1390, Ethanol=590) + >>> bagasse = bst.Stream('bagasse', Water=0.4, Fiber=0.6, total_flow=8e4, units='kg/hr') + >>> with bst.System('sys') as sys: + ... D1 = bst.BinaryDistillation('D1', ins=dilute_ethanol, Lr=0.999, Hr=0.89, k=1.25, LHK=('Ethanol', 'Water')) + ... BT = bst.BoilerTurbogenerator('BT') + ... BT.ins[0] = bagasse + >>> sys.simulate() + >>> BT.show() + BoilerTurbogenerator: BT + ins... + [0] bagasse + phase: 'l', T: 298.15 K, P: 101325 Pa + flow (kmol/hr): Water 1.78e+03 + Cellulose 139 + Hemicellulose 101 + Lignin 79.5 + [1] - + phase: 'l', T: 298.15 K, P: 101325 Pa + flow: 0 + [2] - + phase: 'l', T: 298.15 K, P: 101325 Pa + flow (kmol/hr): Water 488 + [3] - + phase: 'g', T: 288.71 K, P: 101560 Pa + flow: 0 + [4] - + phase: 'l', T: 298.15 K, P: 101325 Pa + flow: 0 + [5] - + phase: 'l', T: 298.15 K, P: 101325 Pa + flow (kmol/hr): Ash 0.567 + [6] - + phase: 'g', T: 298.15 K, P: 101325 Pa + flow (kmol/hr): O2 9.85e+03 + N2 4.23e+04 + outs... + [0] emissions + phase: 'g', T: 394.15 K, P: 101325 Pa + flow (kmol/hr): Water 3.19e+03 + CO2 1.98e+03 + O2 7.84e+03 + N2 4.23e+04 + [1] rejected_water_and_blowdown + phase: 'l', T: 298.15 K, P: 101325 Pa + flow (kmol/hr): Water 488 + [2] ash_disposal + phase: 'l', T: 298.15 K, P: 101325 Pa + flow (kmol/hr): Water 0.00944 + Ash 0.567 + + >>> BT.results() # Steam and electricity are produced, so costs are negative + Boiler turbogenerator Units BT + Electricity Power kW -1.31e+05 + Cost USD/hr -1.02e+04 + Low pressure steam Duty kJ/hr -7.32e+07 + Flow kmol/hr -1.89e+03 + Cost USD/hr -450 + Cooling water Duty kJ/hr -8.42e+07 + Flow kmol/hr 5.75e+04 + Cost USD/hr 28.1 + Fuel (inlet) Flow kg/hr 0 + Cost USD/hr 0 + Ash disposal (outlet) Flow kg/hr 0.737 + Cost USD/hr 0.0234 + Design Work kW 1.33e+05 + Flow rate kg/hr 2.93e+05 + Ash disposal kg/hr 0.737 + Purchase cost Baghouse bags USD 81.1 + Boiler USD 3.33e+07 + Deaerator USD 3.58e+05 + Amine addition pkg USD 4.69e+04 + Hot process water softener system USD 9.16e+04 + Turbogenerator USD 1.94e+07 + Total purchase cost USD 5.32e+07 + Utility cost USD/hr -1.07e+04 + + Notes + ----- + The flow rate of natural gas, lime, and boiler chemicals (streams 3-5) + is set by the BoilerTurbogenerator object during simulation. + + References + ---------- + .. [1] Humbird, D., Davis, R., Tao, L., Kinchin, C., Hsu, D., Aden, A., + Dudgeon, D. (2011). Process Design and Economics for Biochemical + Conversion of Lignocellulosic Biomass to Ethanol: Dilute-Acid + Pretreatment and Enzymatic Hydrolysis of Corn Stover + (No. NREL/TP-5100-47764, 1013269). https://doi.org/10.2172/1013269 + + """ + ticket_name = 'BT' + network_priority = 1 + boiler_blowdown = 0.03 + # Reverse osmosis (RO) typically rejects 25% of water, but the boiler-feed water is assumed to come after RO. + # Setting this parameter to a fraction more than zero effectively assumes that this unit accounts for reverse osmosis. + RO_rejection = 0 + _N_ins = 7 + _N_outs = 3 + _units = {'Flow rate': 'kg/hr', + 'Work': 'kW', + 'Ash disposal': 'kg/hr'} + + def __init__(self, ID='', ins=None, + outs=('emissions', + 'rejected_water_and_blowdown', + 'ash_disposal'), + thermo=None, *, + boiler_efficiency=None, + turbogenerator_efficiency=None, + side_steam=None, + agent=None, + other_agents=None, + fuel_price=None, + natural_gas_price=None, + ash_disposal_price=None, + T_emissions=None, + CO2_emissions_concentration=None, + satisfy_system_electricity_demand=None, + boiler_efficiency_basis=None, + fuel_source=None, + oxygen_rich_gas_composition=None, + ): + if boiler_efficiency_basis is None: boiler_efficiency_basis = 'LHV' + if boiler_efficiency is None: boiler_efficiency = 0.80 + if turbogenerator_efficiency is None: turbogenerator_efficiency = 0.85 + if satisfy_system_electricity_demand is None: satisfy_system_electricity_demand = True + if fuel_source is None: fuel_source = 'CH4' + if oxygen_rich_gas_composition is None: oxygen_rich_gas_composition = dict(O2=21, N2=79, phase='g', units='kg/hr') + if CO2_emissions_concentration is None: CO2_emissions_concentration = 0.055 # Usually between 4 - 7 for biomass and natural gas (https://www.sciencedirect.com/science/article/pii/S0957582021005127) + bst.Facility.__init__(self, ID, ins, outs, thermo) + settings = bst.settings + self.boiler_efficiency_basis = boiler_efficiency_basis + self.agent = agent = agent or settings.get_heating_agent('low_pressure_steam') + self.fuel_source = fuel_source + self.define_utility('Fuel', self.fuel) + self.define_utility('Ash disposal', self.ash_disposal) + self.boiler_efficiency = boiler_efficiency + self.turbogenerator_efficiency = turbogenerator_efficiency + self.steam_utilities = [] + self.steam_demand = agent.to_stream() + self.side_steam = side_steam + self.other_agents = [i for i in settings.heating_agents if i is not agent] if other_agents is None else other_agents + self.CO2_emissions_concentration = CO2_emissions_concentration + self.oxygen_rich_gas_composition = oxygen_rich_gas_composition + if T_emissions is not None: + raise ValueError('setting T_emissions is not yet implemented') + # T_emissions should dictate the efficiency of either the boiler or the turbogenerator. + # Note that T_emissions is tied to the energy balance. + self.T_emissions = T_emissions + if natural_gas_price is not None: + self.fuel_price = natural_gas_price + elif fuel_price is not None: + self.fuel_price = fuel_price + if ash_disposal_price is not None: self.ash_disposal_price = ash_disposal_price + self.satisfy_system_electricity_demand = satisfy_system_electricity_demand + + def _get_desulfurization_rxn_and_coreactant(self): + try: + return self.desulfurization_reaction, self._ID_lime + except: + chemicals = self.chemicals + CAS_lime = '1305-62-0' + has_common_name = 'Ca(OH)2' in chemicals + if CAS_lime in chemicals or has_common_name: + if not has_common_name: chemicals.set_synonym(CAS_lime, 'Ca(OH)2') + self.desulfurization_reaction = rxn = tmo.Reaction( + 'SO2 + Ca(OH)2 + 0.5 O2 -> CaSO4 + H2O', 'SO2', 0.92, chemicals + ) + self._ID_lime = ID = 'Ca(OH)2' + return rxn, ID + CAS_lime = '1305-78-8' + has_common_name = 'CaO' in chemicals + if CAS_lime in chemicals or has_common_name: + if not has_common_name: chemicals.set_synonym(CAS_lime, 'CaO') + self.desulfurization_reaction = rxn = tmo.Reaction( + 'SO2 + CaO + 0.5 O2 -> CaSO4', 'SO2', 0.92, chemicals + ) + self._ID_lime = ID = 'CaO' + return rxn, ID + raise RuntimeError( + "lime is required for boiler, but no chemical 'CaO' or 'Ca(OH)2' " + "available in thermodynamic property package" + ) + + @property + def emissions(self): + return self.outs[0] + + @property + def blowdown_water(self): + return self.outs[1] + + @property + def makeup_water(self): + """[Stream] Makeup water due to boiler blowdown.""" + return self.ins[2] + + @property + def fuel(self): + """[Stream] Fuel used to satisfy steam and electricity requirements.""" + return self.ins[3] + natural_gas = fuel + + @property + def air(self): + """[Stream] Air or oxygen rich gas used to supply oxygen for combustion.""" + return self.ins[6] + oxygen_rich_gas = air + + @property + def ash_disposal(self): + """[Stream] Ash disposal.""" + return self.outs[2] + + @property + def fuel_price(self): + """[Float] Price of natural gas, same as `bst.stream_utility_prices['Natural gas']`.""" + return bst.stream_utility_prices['Fuel'] + @fuel_price.setter + def fuel_price(self, new_price): + bst.stream_utility_prices['Fuel'] = new_price + natural_gas_price = fuel_price + + @property + def ash_disposal_price(self): + """[Float] Price of ash disposal, same as `bst.stream_utility_prices['Ash disposal']`.""" + return bst.stream_utility_prices['Ash disposal'] + + @ash_disposal_price.setter + def ash_disposal_price(self, ash_disposal_price): + bst.stream_utility_prices['Ash disposal'] = ash_disposal_price + + def _run(self): pass + + def _load_utility_agents(self): + steam_utilities = self.steam_utilities + steam_utilities.clear() + agent = self.agent + units = self.other_units + for agent in (*self.other_agents, agent): + ID = agent.ID + for u in units: + for hu in u.heat_utilities: + agent = hu.agent + if agent and agent.ID == ID: + steam_utilities.append(hu) + self.electricity_demand = sum([u.power_utility.consumption for u in units]) + + def _design(self): + B_eff = self.boiler_efficiency + TG_eff = self.turbogenerator_efficiency + steam_demand = self.steam_demand + Design = self.design_results + chemicals = self.chemicals + self._load_utility_agents() + mol_steam = sum([i.flow for i in self.steam_utilities]) + feed_solids, feed_gas, makeup_water, fuel, lime, chems, oxygen_rich_gas = self.ins + oxygen_rich_gas.empty() + if self.fuel_source == 'CH4': + fuel.phase = 'g' + fuel.set_property('T', 60, 'degF') + fuel.set_property('P', 14.73, 'psi') + emissions, blowdown_water, ash_disposal = self.outs + if not lime.price: + lime.price = 0.19937504680689402 + if not chems.price: + chems.price = 4.995862254032183 + H_steam = sum([i.duty for i in self.steam_utilities]) + side_steam = self.side_steam + if side_steam: + H_steam += side_steam.H + mol_steam += side_steam.F_mol + steam_demand.imol['7732-18-5'] = mol_steam + self.combustion_reactions = combustion_rxns = chemicals.get_combustion_reactions() + non_empty_feeds = [i for i in (feed_solids, feed_gas) if not i.isempty()] + boiler_efficiency_basis = self.boiler_efficiency_basis + fuel_source = self.fuel_source + def calculate_excess_electricity_at_natual_gas_flow(fuel_flow): + if fuel_flow: + fuel_flow = abs(fuel_flow) + fuel.imol[fuel_source] = fuel_flow + else: + fuel.empty() + if boiler_efficiency_basis == 'LHV': + H_combustion = fuel.LHV + for feed in non_empty_feeds: H_combustion += feed.LHV + elif boiler_efficiency_basis == 'HHV': + H_combustion = fuel.HHV + for feed in non_empty_feeds: H_combustion += feed.HHV + else: + raise ValueError( + f"invalid boiler efficiency basis {boiler_efficiency_basis}; " + f"valid values include 'LHV', or 'HHV'" + ) + self.H_content = H_content = B_eff * H_combustion + self.H_loss_to_emissions = H_combustion - H_content + H_electricity = H_content - H_steam # Heat available for the turbogenerator + electricity = H_electricity * TG_eff # Electricity produced + self.cooling_duty = electricity - H_electricity + Design['Work'] = work = electricity / 3600 + duty_over_mol = 39000 # kJ / mol-superheated steam + self.total_steam = H_content / duty_over_mol #: [float] Total steam produced by the boiler (kmol/hr) + Design['Flow rate'] = flow_rate = self.total_steam * 18.01528 + + if self.satisfy_system_electricity_demand: + boiler = self.cost_items['Boiler'] + rate_boiler = boiler.kW * flow_rate / boiler.S + return work - self.electricity_demand - rate_boiler + else: + return work + + self._excess_electricity_without_fuel = excess_electricity = calculate_excess_electricity_at_natual_gas_flow(0) + if excess_electricity < 0: + f = calculate_excess_electricity_at_natual_gas_flow + lb = 0. + fuel.imol[fuel_source] = 1 + ub = - excess_electricity * 3600 / fuel.LHV + while f(ub) < 0.: + lb = ub + ub *= 2 + flx.IQ_interpolation(f, lb, ub, xtol=1, ytol=1) + + if self.cooling_duty > 0.: + # In the event that no electricity is produced and the solver + # solution for natural gas is slightly below the requirement for steam + # (this would lead to a positive duty). + self.cooling_duty = 0. + Design['Work'] = 0. + + emissions.T = 298.15 # Will be updated later with the energy balance + emissions.P = 101325 + emissions.phase = 'g' + emissions_mol = emissions.mol + emissions_mol[:] = fuel.mol + for feed in non_empty_feeds: emissions_mol[:] += feed.mol + combustion_rxns.force_reaction(emissions_mol) + O2_consumption = -emissions.imol['O2'] + oxygen_rich_gas.reset_flow(**self.oxygen_rich_gas_composition) + z_O2 = oxygen_rich_gas.imol['O2'] / oxygen_rich_gas.F_mol + oxygen_rich_gas.F_mol = O2_consumption / z_O2 + emissions_mol += oxygen_rich_gas.mol + F_emissions = emissions.F_mass + z_CO2 = emissions.imass['CO2'] / F_emissions + z_CO2_target = self.CO2_emissions_concentration + if z_CO2 > z_CO2_target: + F_emissions_new = z_CO2 * F_emissions / z_CO2_target + dF_emissions = F_emissions_new - F_emissions + oxygen_rich_gas.F_mass = F_mass_O2_new = oxygen_rich_gas.F_mass + dF_emissions + emissions_mol += oxygen_rich_gas.mol * (dF_emissions / F_mass_O2_new) + emissions.H += self.H_loss_to_emissions + hu_cooling = bst.HeatUtility() + hu_cooling(self.cooling_duty, steam_demand.T) + hus_heating = bst.HeatUtility.sum_by_agent(self.steam_utilities) + for hu in hus_heating: hu.reverse() + self.heat_utilities = [*hus_heating, hu_cooling] + water_index = chemicals.index('7732-18-5') + makeup_water.mol[water_index] = blowdown_water.mol[water_index] = ( + self.total_steam * self.boiler_blowdown * 1 / (1 - self.RO_rejection) + ) + ash_IDs = [i.ID for i in self.chemicals if not i.formula] + emissions_mol = emissions.mol + SO2_produced = 0 + if 'SO2' in chemicals: + SO2_produced += emissions.imol['SO2'] + if 'CaSO4' in chemicals: + SO2_produced += emissions.imol['CaSO4'] + ash_IDs.append('CaSO4') + if SO2_produced: + rxn, ID_lime = self._get_desulfurization_rxn_and_coreactant() + # FGD lime scaled based on SO2 generated, + # 20% stoichiometric excess based on P52 of ref [1] + rxn.force_reaction(emissions) + lime.imol[ID_lime] = lime_mol = SO2_produced * 1.2 + emissions_mol.remove_negatives() + else: + lime.empty() + # About 0.4536 kg/hr of boiler chemicals are needed per 234484 kg/hr steam produced + chems.imol['Ash'] = boiler_chems = 1.9345e-06 * Design['Flow rate'] + ash_disposal.empty() + ash_disposal.copy_flow(emissions, IDs=tuple(ash_IDs), remove=True) + ash_disposal.imol['Ash'] += boiler_chems + dry_ash = ash_disposal.F_mass + moisture = min(emissions.imass['Water'], dry_ash * 0.3) # ~20% moisture + ash_disposal.imass['Water'] = moisture + emissions.imass['Water'] -= moisture + Design['Ash disposal'] = dry_ash + moisture + if SO2_produced: + if ID_lime == 'Ca(OH)2': # Ca(OH)2 + lime.imol['Water'] = 4 * lime_mol # Its a slurry + else: # CaO + lime.imol['Water'] = 5 * lime_mol + + def _cost(self): + self._decorated_cost() + self.power_utility.production = self.design_results['Work'] + + +# %% +# ============================================================================= +# Custom Heat Exchanger Network with minimum external heat utility constraint +# ============================================================================= + +class HXNWithMin(bst.facilities.HeatExchangerNetwork): + """Heat-exchanger network that caps maximum heat recovery so that the required + external hot utility duty is not lower than ``min_heat_util`` (kJ/hr). + + Parameters + ---------- + min_heat_util : float, optional + Minimum positive duty (kJ/hr) that must still be met by external + heating utilities (e.g., steam). Defaults to 0 (i.e., no cap). + T_min_app : float, optional + Minimum approach temperature for pinch analysis [K]. + """ + + def __init__(self, ID='', *, T_min_app=10.0, min_heat_util=0.0): + self.min_heat_util = min_heat_util # kJ/hr + super().__init__(ID, T_min_app=T_min_app) + + def _cost(self): + # First run default costing which also computes heat utility loads + super()._cost() + # After synthesis, `actual_heat_util_load` is positive if external hot utility is required + # and negative/zero if process heat is sufficient. + hot_util_load = getattr(self, 'actual_heat_util_load', None) + if hot_util_load is None: + return # Nothing to do + if hot_util_load < self.min_heat_util: + make_up_duty = self.min_heat_util - hot_util_load # kJ/hr to be supplied by steam + # Create additional heat utility demand so external steam is still counted + hu = bst.HeatUtility.get_agent('low_pressure_steam') + hu(make_up_duty, 423.15) # Typical LP steam temperature ~150 °C (423 K) + self.heat_utilities.append(hu) + + diff --git a/biorefineries/microalgae/lca.py b/biorefineries/microalgae/lca.py new file mode 100644 index 00000000..e57410e5 --- /dev/null +++ b/biorefineries/microalgae/lca.py @@ -0,0 +1,446 @@ +# -*- coding: utf-8 -*- +""" +Created on Sat July 20 10:50:00 2025 + +Microalgae biorefinery to produce medium chain fatty acids +by anaerobic fermentation without external electron donor addition- LCA analysis + +References +---------- +[1] BioSTEAM Documentation: + https://biosteam.readthedocs.io/en/latest/tutorial/Creating_a_System.html +[2] Cortes-Peña et al., BioSTEAM: A Fast and Flexible Platform for the Design, + Simulation, and Techno-Economic Analysis of Biorefineries under Uncertainty. + ACS Sustainable Chem. Eng. 2020, 8 (8), 3302–3310. +[3] 3-Hydroxypropionic acid biorefineries project: + https://github.com/BioSTEAMDevelopmentGroup/Bioindustrial-Park/tree/master/biorefineries/HP +[4] Succinic projest + https://github.com/BioSTEAMDevelopmentGroup/Bioindustrial-Park/tree/master/biorefineries/succinic + +@author: Xingdong Shi +@version: 0.0.9 +""" + +from thermosteam import Stream +import numpy as np +import pandas as pd +import biosteam as bst +from .system import create_microalgae_MCCA_production_sys +from ._chemicals import chems + +__all__ = ['LCA'] + +class LCA: + GWP_key = 'GWP_100' + FEC_key = 'FEC' + + def __init__(self, + system, + CFs, + main_product, + main_product_chemical_IDs, + boiler, + by_products=[], + complex_feeds={}, + cooling_tower=None, + chilled_water_processing_units=[], + functional_unit=None, + functional_quantity_per_h_fn=None, + add_EOL_GWP=False, + input_biogenic_carbon_streams=[]): + + self.system = system + system._LCA = self + self.chemicals = chemicals = main_product.chemicals + self.units = self.system.units + self.flowsheet = self.system.flowsheet + self.streams = self.system.streams + self.input_biogenic_carbon_streams = input_biogenic_carbon_streams + self.add_EOL_GWP = add_EOL_GWP + self.complex_feeds = complex_feeds + + kg_product_per_h_fn = lambda: self.main_product.imass[self.main_product_chemical_IDs].sum() + MJ_LHV_product_per_h_fn = lambda: self.main_product.LHV + MJ_HHV_product_per_h_fn = lambda: self.main_product.HHV + if not functional_unit: + self.functional_unit = functional_unit = '1 kg' + self.functional_quantity_per_h_fn = kg_product_per_h_fn + + elif not functional_quantity_per_h_fn: + self.functional_unit = functional_unit + if functional_unit == '1 kg': + self.functional_quantity_per_h_fn = kg_product_per_h_fn + elif functional_unit == '1 metric ton': + self.functional_quantity_per_h_fn = lambda: kg_product_per_h_fn() * 1000. + elif functional_unit == '1 MJ (by LHV)': + self.functional_quantity_per_h_fn = MJ_LHV_product_per_h_fn + elif functional_unit == '1 MJ (by HHV)': + self.functional_quantity_per_h_fn = MJ_HHV_product_per_h_fn + + self.CFs = CFs + self._chemical_IDs = [chem.ID for chem in chemicals] + self._CF_streams = _CF_streams = {} + for impact_category in CFs.keys(): + _CF_streams[impact_category] = ic_CF_stream = Stream(f'{system.ID}_{impact_category}_CF_stream') + for k, v in CFs[impact_category].items(): + if k not in complex_feeds and k != 'Electricity': + try: + ic_CF_stream.imass[k] = v + except: + pass + + self._LCA_stream = Stream(f'{system.ID}_LCA_stream') + self.main_product_chemical_IDs = main_product_chemical_IDs + self.main_product = main_product + self.by_products = by_products + self.chem_IDs = [i.ID for i in chemicals] + self.CT = self.cooling_tower = cooling_tower + self.CWP_units = self.chilled_water_processing_units = chilled_water_processing_units + self._CO2_MW = self.chemicals.CO2.MW + + @property + def products(self): + return [self.main_product] + self.by_products + + @property + def functional_quantity_per_h(self): + return self.functional_quantity_per_h_fn() + + @property + def emissions(self): + emissions = list(self.system.products) + for i in self.products: + if i in emissions: + emissions.remove(i) + return emissions + + @property + def feeds(self): + return self.system.feeds + + @property + def LCA_stream(self): + _LCA_stream = self._LCA_stream + to_mix = list(self.feeds) + for s, m_k in self.complex_feeds.values(): + if s in to_mix: + to_mix.remove(s) + _LCA_stream.mix_from(to_mix) + return _LCA_stream + + def get_material_impact_array(self, impact_category): + return self.LCA_stream.mass*self._CF_streams[impact_category].mass + + def get_material_impact(self, impact_category): + return self.get_material_impact_array(impact_category).sum() / self.functional_quantity_per_h + + @property + def net_electricity(self): + return self.system.power_utility.rate + + def get_net_electricity_impact(self, impact_category): + if 'Electricity' in self.CFs[impact_category]: + return self.net_electricity * self.CFs[impact_category]['Electricity'] / self.functional_quantity_per_h + return 0 + + @property + def EOL_GWP(self): + return sum([i.get_atomic_flow('C') for i in [self.main_product] + self.by_products]) *\ + self._CO2_MW/self.functional_quantity_per_h + + @property + def direct_emissions_GWP(self): + return sum([stream.get_atomic_flow('C') for stream in self.emissions]) *\ + self._CO2_MW / self.functional_quantity_per_h + + @property + def biogenic_emissions_GWP(self): + return sum([i.get_atomic_flow('C') for i in self.input_biogenic_carbon_streams]) *\ + self._CO2_MW/self.functional_quantity_per_h + + @property + def direct_non_biogenic_emissions_GWP(self): + return self.direct_emissions_GWP -\ + self.biogenic_emissions_GWP +\ + int(self.add_EOL_GWP)*self.EOL_GWP + + def get_complex_feeds_impact(self, impact_category): + tot_cfs_impact = 0. + impact_CFs = self.CFs[impact_category] + for k, (s, mass_kind) in self.complex_feeds.items(): + if k in impact_CFs: + mass = s.F_mass + if mass_kind=='dry': mass -= s.imass['Water'] + tot_cfs_impact += impact_CFs[k] * mass + return tot_cfs_impact/self.functional_quantity_per_h + + def get_total_impact(self, impact_category): + tot_impact = self.get_complex_feeds_impact(impact_category) +\ + self.get_material_impact(impact_category) +\ + self.get_net_electricity_impact(impact_category) + if impact_category in ('GWP', 'GWP_100', 'GWP100'): + tot_impact += self.direct_non_biogenic_emissions_GWP + return tot_impact + + + def get_material_impact_breakdown(self, impact_category): + chemical_impact_dict = {'H2SO4':0, 'NaOH':0, 'NH4OH':0, 'CH4':0, 'CO2':0, 'OleylAlcohol':0} + LCA_stream = self.LCA_stream + ic_CF_stream = self._CF_streams[impact_category] + functional_quantity_per_h = self.functional_quantity_per_h + chemical_impact_dict_additional = {ID: LCA_stream.imass[ID] * ic_CF_stream.imass[ID] / functional_quantity_per_h + for ID in self.chem_IDs} + for k in list(chemical_impact_dict_additional.keys()): + if chemical_impact_dict_additional[k] == 0.: del(chemical_impact_dict_additional[k]) + chemical_impact_dict.update(chemical_impact_dict_additional) + return chemical_impact_dict + + @property + def GWP(self): + return self.get_total_impact(self.GWP_key) + + @property + def material_FEC(self): + return self.get_material_impact(self.FEC_key) + + @property + def feedstock_FEC(self): + return self.get_complex_feeds_impact(self.FEC_key) + + @property + def net_electricity_FEC(self): + return self.get_net_electricity_impact(self.FEC_key) + + @property + def material_FEC_breakdown(self): + return self.get_material_impact_breakdown(self.FEC_key) + + @property + def FEC(self): + return self.get_total_impact(self.FEC_key) + + def get_detailed_GWP_breakdown(self): + breakdown = {} + gwp_key = self.GWP_key + gwp_cfs = self.CFs[gwp_key] + functional_quantity = self.functional_quantity_per_h + + # 1. Feedstock impact + feedstock_impact = self.get_complex_feeds_impact(gwp_key) + breakdown['feedstock'] = feedstock_impact + + # 2. Material impacts from various chemicals + chemical_impacts = self.get_material_impact_breakdown(gwp_key) + + # Extract major chemicals (enzymes and oleyl alcohol go to other materials) + major_chemicals = ['NaOH', 'NH4OH', 'Lime'] + + for chem in major_chemicals: + if chem in chemical_impacts and abs(chemical_impacts[chem]) > 1e-6: + breakdown[chem.lower()] = chemical_impacts[chem] + + # 3. Other materials (including enzymes and oleyl alcohol) + other_materials = 0.0 + for chem_id, impact in chemical_impacts.items(): + if chem_id not in major_chemicals and abs(impact) > 1e-6: + other_materials += impact + breakdown['other_materials'] = other_materials + + + + # 5. Net electricity impact + net_electricity_impact = self.get_net_electricity_impact(gwp_key) + breakdown['net_electricity'] = net_electricity_impact + + # 6. Direct emissions - only show final non-biogenic emissions + breakdown['direct_non_biogenic_emissions'] = self.direct_non_biogenic_emissions_GWP + + breakdown['total_actual'] = self.GWP + + return breakdown + + def print_GWP_breakdown(self, sort_by_magnitude=True, show_small_contributions=False): + breakdown = self.get_detailed_GWP_breakdown() + + print(f"Total GWP: {breakdown['total_actual']:.4f} kg CO2e/{self.functional_unit}") + print("\nComponent contributions:") + print("-" * 80) + print(f"{'Component Name':<30} {'Contribution':<15} {'Percentage':<10} {'Unit'}") + print("-" * 80) + + display_items = [] + total_gwp = breakdown['total_actual'] + component_map = { + 'feedstock': 'Microalgae feedstock', + 'naoh': 'Sodium hydroxide (NaOH)', + 'nh4oh': 'Ammonium hydroxide (NH4OH)', + 'lime': 'Lime', + 'other_materials': 'Other materials', + 'net_electricity': 'Net electricity', + 'direct_non_biogenic_emissions': 'Direct non-biogenic emissions' + } + + for key, value in breakdown.items(): + if key == 'total_actual': + continue + + if not show_small_contributions and abs(value) < 1e-4: + continue + + display_name = component_map.get(key, key) + percentage = (value / total_gwp * 100) if total_gwp != 0 else 0 + display_items.append((display_name, value, percentage)) + + if sort_by_magnitude: + display_items.sort(key=lambda x: abs(x[1]), reverse=True) + + for name, value, percentage in display_items: + unit = f"kg CO2e/{self.functional_unit}" + print(f"{name:<30} {value:>14.4f} {percentage:>9.1f}% {unit}") + + return breakdown + + def get_GWP_breakdown_dataframe(self): + breakdown = self.get_detailed_GWP_breakdown() + + component_map = { + 'feedstock': 'Microalgae feedstock', + 'naoh': 'Sodium hydroxide (NaOH)', + 'nh4oh': 'Ammonium hydroxide (NH4OH)', + 'lime': 'Lime', + 'other_materials': 'Other materials', + 'net_electricity': 'Net electricity', + 'direct_non_biogenic_emissions': 'Direct non-biogenic emissions' + } + + data = [] + total_gwp = breakdown['total_actual'] + + for key, value in breakdown.items(): + if key == 'total_actual': + continue + + display_name = component_map.get(key, key) + percentage = (value / total_gwp * 100) if total_gwp != 0 else 0 + + data.append({ + 'Component': display_name, + 'Original_key': key, + f'GWP_contribution (kg CO2e/{self.functional_unit})': value, + 'Percentage (%)': percentage, + 'Absolute_value': abs(value) + }) + + df = pd.DataFrame(data) + df = df.sort_values('Absolute_value', ascending=False) + df = df.drop('Absolute_value', axis=1) + + return df + + +def create_microalgae_lca(system, main_product, main_product_chemical_IDs, boiler): + CFs = { + # ============================================================================= + # 100-year global warming potential (GWP) in kg CO2-eq/kg + # ============================================================================= + 'GWP_100':{ + 'Electricity': 0.36, # [kg*CO2*eq / kWhr] From GREET; NG-Fired Simple-Cycle Gas Turbine CHP Plant + # 0.66 is the GWP from producing diesel from GREET; Conventional diesel from crude oil for US Refineries. + # Downstream fuel emissions are added in. Accounts for how biodiesel has less energy than diesel. + 'CH4': 0.33, # Natural gas from shell conventional recovery, GREET; includes non-biogenic emissions + 'H2SO4': 0.04447, # kg CO2e/kg + 'NaOH': 2.01, # GREET + 'NH4OH': 1.28304, # multiplied by chemicals.NH3.MW/chemicals.NH4OH.MW, + 'CalciumDihydroxide': 1.29, # GREET + 'Ethanol': 1.44, # from BDO project + + 'GlucoAmylase': 6.16, # from Succinic project + 'AlphaAmylase': 6.16, # from Succinic project + + 'OleylAlcohol': 1.9, # kg CO2e/kg + 'Lime': 1.29 * 56.0774/74.093, # CaO to Ca(OH)2 + 'Microalgae': 0.1, # + 'CO2': 0.87104, # ecoinvent 3.8 carbon dioxide production, liquid, RoW + 'CH4': 0.40, # NA NG from shale and conventional recovery + 'Electricity': 0.4490 # kg CO2-eq/kWh GREET 2022 US Mix # assume production==consumption, both in kg CO2-eq/kWh + }, + + # ============================================================================= + # Fossil energy consumption (FEC), in MJ/kg of material + # ============================================================================= + 'FEC': { + 'Electricity': 10.0, # MJ/kWh + 'H2SO4': 568.98/1e3, + 'NaOH': 29, # MJ/kg + 'NH4OH': 42 * 0.4860, # multiplied by chemicals.NH3.MW/chemicals.NH4OH.MW, + 'CalciumDihydroxide': 4.0, # MJ/kg + 'CH4': 55.0, # MJ/kg + 'Ethanol': 16, + 'OleylAlcohol': 45.0, # MJ/kg + 'GlucoAmylase': 90.0, # MJ/kg + 'AlphaAmylase': 90.0, # MJ/kg + 'Lime': 4.896 * 56.0774/74.093, # CaO to Ca(OH)2 + 'Microalgae': 2.0, # MJ/kg + 'Electricity': 5.724, # MJ/kWh # GREET 2022 US Mix #assume production==consumption, both in MJ/kWh + 'CH4': 50, # NA NG from shale and conventional recovery + } + } + + cooling_tower = next((unit for unit in system.units + if unit.__class__.__name__ == 'CoolingTower'), None) + chilled_water_processing_units = [unit for unit in system.units + if unit.__class__.__name__ == 'ChilledWaterPackage'] + boiler = system.flowsheet.unit.BT601 + + complex_feeds = { + 'Microalgae': (system.feeds[0], 'dry'), + } + + microalgae_feed = next((feed for feed in system.feeds + if 'microalgae' in feed.ID.lower()), + next((feed for feed in system.feeds + if feed.get_atomic_flow('C') > 0), None)) + + input_biogenic_carbon_streams = [microalgae_feed] if microalgae_feed else [] + microalgae_lca = LCA( + system=system, + CFs=CFs, + main_product=main_product, + main_product_chemical_IDs=main_product_chemical_IDs, + boiler=boiler, + by_products=[system.flowsheet.stream.butyric_acid_product, + system.flowsheet.stream.heptanoic_acid_product, + system.flowsheet.stream.caprylic_acid_product], + complex_feeds=complex_feeds, + cooling_tower=cooling_tower, + chilled_water_processing_units=chilled_water_processing_units, + functional_unit='1 kg', + add_EOL_GWP=True, + input_biogenic_carbon_streams=input_biogenic_carbon_streams + ) + + return microalgae_lca + +if __name__ == "__main__": + bst.settings.set_thermo(chems) + microalgae_sys = create_microalgae_MCCA_production_sys() + microalgae_sys.simulate() + main_product = microalgae_sys.flowsheet.stream.caproic_acid_product + main_product_chemical_IDs = ['CaproicAcid'] + boiler = microalgae_sys.flowsheet.unit.BT601 + lca = create_microalgae_lca( + system=microalgae_sys, + main_product=main_product, + main_product_chemical_IDs=main_product_chemical_IDs, + boiler=boiler + ) + + # Call detailed breakdown method + breakdown_data = lca.print_GWP_breakdown(sort_by_magnitude=True, show_small_contributions=False) + # Save DataFrame to CSV file (optional) + # breakdown_df.to_csv('microalgae_gwp_breakdown.csv', index=False, encoding='utf-8-sig') + + + + + diff --git a/biorefineries/microalgae/lca_output.py b/biorefineries/microalgae/lca_output.py new file mode 100644 index 00000000..2208c1df --- /dev/null +++ b/biorefineries/microalgae/lca_output.py @@ -0,0 +1,380 @@ +# -*- coding: utf-8 -*- +""" +Microalgae LCA Output Module + +基于succinic项目的lca_output.py模式，为microalgae项目提供LCA分析、报告和数据导出功能。 + +@author: Xingdong Shi +@version: 0.0.1 +""" + +import sys +import os +import numpy as np +import pandas as pd +import matplotlib.pyplot as plt +from pathlib import Path + +# 添加项目路径以便导入 +current_dir = Path(__file__).parent +project_root = current_dir.parent +sys.path.insert(0, str(project_root)) + +try: + from . import system + from ._chemicals import chems + from .lca import LCA, create_microalgae_lca +except ImportError: + # 如果相对导入失败，尝试绝对导入 + try: + from microalgae import system + from microalgae._chemicals import chems + from microalgae.lca import LCA, create_microalgae_lca + except ImportError: + print("警告：无法导入microalgae模块，请确保在正确的环境中运行") + +def get_gwp_breakdown_data(lca): + """ + 获取GWP分解数据，按照microalgae项目的计算方式 + """ + material_GWP_breakdown = lca.get_material_impact_breakdown('GWP_100') + + gwp_breakdown = { + 'feedstock*': lca.feedstock_GWP, + 'H2SO4': material_GWP_breakdown.get('H2SO4', 0), + 'NaOH': material_GWP_breakdown.get('NaOH', 0), + 'NH4OH': material_GWP_breakdown.get('NH4OH', 0), + 'CalciumDihydroxide': material_GWP_breakdown.get('CalciumDihydroxide', 0), + 'Ethanol': material_GWP_breakdown.get('Ethanol', 0), + 'Octanol': material_GWP_breakdown.get('Octanol', 0), + 'GlucoAmylase': material_GWP_breakdown.get('GlucoAmylase', 0), + 'AlphaAmylase': material_GWP_breakdown.get('AlphaAmylase', 0), + 'CH4': material_GWP_breakdown.get('CH4', 0), + 'net electricity': lca.net_electricity_GWP, + 'direct non-biogenic emissions': lca.direct_non_biogenic_emissions_GWP, + } + + # 计算正值的总和，然后归一化 + tot_positive_GWP = sum([v for v in gwp_breakdown.values() if v > 0]) + if tot_positive_GWP > 0: + for k, v in gwp_breakdown.items(): + gwp_breakdown[k] = v / tot_positive_GWP + else: # Fallback if no positive values + tot_abs_GWP = sum([abs(v) for v in gwp_breakdown.values()]) + if tot_abs_GWP > 0: + for k, v in gwp_breakdown.items(): + gwp_breakdown[k] = v / tot_abs_GWP + + return gwp_breakdown + +def export_gwp_breakdown_data(lca, output_file=None): + """ + 导出GWP分解数据到CSV文件 + + 参数 + ---------- + lca : LCA + LCA对象 + output_file : str, 可选 + 输出文件路径，默认为'microalgae_gwp_breakdown.csv' + """ + if output_file is None: + output_file = 'microalgae_gwp_breakdown.csv' + + gwp_breakdown = get_gwp_breakdown_data(lca) + + # 创建DataFrame + df = pd.DataFrame(list(gwp_breakdown.items()), columns=['Component', 'Percentage']) + df['Percentage'] = df['Percentage'] * 100 # 转换为百分比 + + # 保存到CSV + df.to_csv(output_file, index=False) + print(f"GWP分解数据已保存到: {output_file}") + + return df + +def print_microalgae_lca_results(lca): + """ + 打印microalgae LCA结果 + + 参数 + ---------- + lca : LCA + LCA对象 + """ + print("\n" + "="*60) + print("MICROALGAE 生物炼制厂生命周期评估结果") + print("="*60) + + print(f"\n功能单位: {lca.functional_unit}") + print(f"系统碳平衡: {lca.system_carbon_balance:.4f}") + + # 主要影响类别结果 + print(f"\n----- 主要影响类别 -----") + print(f"GWP (kg CO2e/{lca.functional_unit}): {lca.GWP:.4f}") + print(f"FEC (MJ/{lca.functional_unit}): {lca.FEC:.4f}") + print(f"WC (m³/{lca.functional_unit}): {lca.WC:.4f}") + + # GWP分解 + print(f"\n----- GWP分解 -----") + print(f"物质影响: {lca.material_GWP:.4f} ({lca.material_GWP/lca.GWP*100 if lca.GWP else 0:.1f}%)") + print(f"原料影响: {lca.feedstock_GWP:.4f} ({lca.feedstock_GWP/lca.GWP*100 if lca.GWP else 0:.1f}%)") + print(f"电力影响: {lca.net_electricity_GWP:.4f} ({lca.net_electricity_GWP/lca.GWP*100 if lca.GWP else 0:.1f}%)") + print(f"直接排放(非生物源): {lca.direct_non_biogenic_emissions_GWP:.4f} ({lca.direct_non_biogenic_emissions_GWP/lca.GWP*100 if lca.GWP else 0:.1f}%)") + + # 物质影响明细 + print(f"\n----- 物质影响明细 (kg CO2e/{lca.functional_unit}) -----") + material_breakdown = lca.get_material_impact_breakdown('GWP_100') + for component, impact in material_breakdown.items(): + if abs(impact) > 1e-6: + print(f"{component:<20}: {impact:8.4f}") + + # 电力影响详情 + print(f"\n----- 电力影响详情 -----") + elec_CF = lca.CFs['GWP_100']['Electricity'] + net_kWh = lca.net_electricity + per_unit = net_kWh * elec_CF / lca.functional_quantity_per_h + print(f"净电力消耗: {net_kWh:.3f} kWh/h") + print(f"电力影响因子: {elec_CF} kg CO2e/kWh") + print(f"电力影响: {per_unit:.4f} kg CO2e/{lca.functional_unit}") + + # 原料影响详情 + print(f"\n----- 原料影响详情 -----") + for ID, (stream, mass_kind) in lca.complex_feeds.items(): + impact = lca.get_complex_feed_impact_by_ID('GWP_100', ID) + print(f"{ID:<20}: {impact:8.4f} kg CO2e/{lca.functional_unit}") + + # 直接排放详情 + print(f"\n----- 直接排放详情 -----") + print(f"总直接排放GWP: {lca.direct_emissions_GWP:.4f}") + print(f"生物来源排放GWP: {lca.biogenic_emissions_GWP:.4f}") + print(f"非生物来源直接排放GWP: {lca.direct_non_biogenic_emissions_GWP:.4f}") + + # 排放流详细 + print(f"\n排放流详细:") + for stream in lca.emissions: + C = stream.get_atomic_flow("C") + if C > 0: + stream_gwp = C * 44 / lca.functional_quantity_per_h + if abs(stream_gwp) > 1e-6: + print(f" {stream.ID:<20}: {stream_gwp:8.4f} kg CO2e/{lca.functional_unit}") + + # GWP分解数据 + print(f"\n----- GWP分解数据 -----") + gwp_breakdown = get_gwp_breakdown_data(lca) + + print("GWP分解数据:") + print("="*50) + total_percentage = 0 + for component, percentage in gwp_breakdown.items(): + print(f"{component:<30} {percentage*100:6.2f}%") + if percentage > 0: # 只计算正值 + total_percentage += percentage * 100 + print("-"*50) + print(f"{'总计':<30} {total_percentage:6.2f}%") + + print("\n" + "="*60) + +def create_microalgae_lca_scenario(system, main_product, main_product_chemical_IDs, boiler=None): + """ + 创建microalgae LCA场景 + + 参数 + ---------- + system : System + 微藻生产系统 + main_product : Stream + 主产品流 + main_product_chemical_IDs : list + 主产品化学品ID列表 + boiler : Unit, 可选 + 锅炉设备单元 + + 返回 + ---------- + LCA + 生命周期评估对象 + """ + # 创建LCA对象 + lca = create_microalgae_lca( + system=system, + main_product=main_product, + main_product_chemical_IDs=main_product_chemical_IDs, + boiler=boiler, + add_EOL_GWP=True + ) + + return lca + +def export_lca_results_to_excel(lca, output_file=None): + """ + 导出LCA结果到Excel文件 + + 参数 + ---------- + lca : LCA + LCA对象 + output_file : str, 可选 + 输出文件路径，默认为'microalgae_lca_results.xlsx' + """ + if output_file is None: + output_file = 'microalgae_lca_results.xlsx' + + # 创建结果字典 + results = { + '基本信息': { + '功能单位': [lca.functional_unit], + '系统碳平衡': [lca.system_carbon_balance], + '主产品': [lca.main_product.ID], + }, + '主要影响类别': { + 'GWP (kg CO2e/功能单位)': [lca.GWP], + 'FEC (MJ/功能单位)': [lca.FEC], + 'WC (m³/功能单位)': [lca.WC], + }, + 'GWP分解': { + '物质影响 (kg CO2e/功能单位)': [lca.material_GWP], + '原料影响 (kg CO2e/功能单位)': [lca.feedstock_GWP], + '电力影响 (kg CO2e/功能单位)': [lca.net_electricity_GWP], + '直接排放(非生物源) (kg CO2e/功能单位)': [lca.direct_non_biogenic_emissions_GWP], + }, + '物质影响明细': lca.get_material_impact_breakdown('GWP_100'), + 'GWP分解百分比': get_gwp_breakdown_data(lca) + } + + # 创建Excel文件 + with pd.ExcelWriter(output_file, engine='openpyxl') as writer: + # 基本信息 + pd.DataFrame(results['基本信息']).to_excel(writer, sheet_name='基本信息', index=False) + + # 主要影响类别 + pd.DataFrame(results['主要影响类别']).to_excel(writer, sheet_name='主要影响类别', index=False) + + # GWP分解 + pd.DataFrame(results['GWP分解']).to_excel(writer, sheet_name='GWP分解', index=False) + + # 物质影响明细 + material_df = pd.DataFrame(list(results['物质影响明细'].items()), + columns=['物质', '影响 (kg CO2e/功能单位)']) + material_df.to_excel(writer, sheet_name='物质影响明细', index=False) + + # GWP分解百分比 + gwp_percent_df = pd.DataFrame(list(results['GWP分解百分比'].items()), + columns=['组分', '百分比']) + gwp_percent_df['百分比'] = gwp_percent_df['百分比'] * 100 + gwp_percent_df.to_excel(writer, sheet_name='GWP分解百分比', index=False) + + print(f"LCA结果已导出到: {output_file}") + return output_file + +def compare_lca_scenarios(scenarios, output_file=None): + """ + 比较多个LCA场景 + + 参数 + ---------- + scenarios : dict + 场景字典，格式: {scenario_name: lca_object} + output_file : str, 可选 + 输出文件路径 + """ + if output_file is None: + output_file = 'microalgae_lca_comparison.xlsx' + + # 创建比较数据 + comparison_data = {} + + for scenario_name, lca in scenarios.items(): + comparison_data[scenario_name] = { + 'GWP (kg CO2e/功能单位)': lca.GWP, + 'FEC (MJ/功能单位)': lca.FEC, + 'WC (m³/功能单位)': lca.WC, + '物质影响 (kg CO2e/功能单位)': lca.material_GWP, + '原料影响 (kg CO2e/功能单位)': lca.feedstock_GWP, + '电力影响 (kg CO2e/功能单位)': lca.net_electricity_GWP, + '直接排放(非生物源) (kg CO2e/功能单位)': lca.direct_non_biogenic_emissions_GWP, + } + + # 创建DataFrame并保存 + df = pd.DataFrame(comparison_data).T + df.to_excel(output_file) + + print(f"场景比较结果已保存到: {output_file}") + return df + +def generate_lca_report(lca, output_file=None): + """ + 生成详细的LCA报告 + + 参数 + ---------- + lca : LCA + LCA对象 + output_file : str, 可选 + 输出文件路径 + """ + if output_file is None: + output_file = 'microalgae_lca_report.txt' + + report = lca.generate_report(output_file) + + if output_file: + print(f"详细报告已保存到: {output_file}") + + return report + +# 主程序 +if __name__ == "__main__": + import biosteam as bst + + # 初始化系统 + bst.settings.set_thermo(chems) + microalgae_sys = system.create_microalgae_MCCA_production_sys() + microalgae_sys.simulate() + + # 获取主产品流和主产品化学品ID + main_product = microalgae_sys.flowsheet.stream.caproic_acid_product + main_product_chemical_IDs = ['CaproicAcid'] + + # 尝试找到锅炉涡轮发电机 + boiler = None + for unit in microalgae_sys.units: + # 优先查找锅炉涡轮发电机 + if 'BT' in unit.ID and hasattr(unit, 'turbogenerator_efficiency'): + boiler = unit + break + + # 如果没有找到锅炉涡轮发电机，查找普通锅炉 + if boiler is None: + for unit in microalgae_sys.units: + if hasattr(unit, 'heat_utilities') and unit.heat_utilities: + boiler = unit + break + + # 如果没有找到锅炉，使用第一个热交换器 + if boiler is None: + for unit in microalgae_sys.units: + if unit.__class__.__name__ == 'HXutility': + boiler = unit + break + + # 创建LCA对象 + lca = create_microalgae_lca_scenario( + system=microalgae_sys, + main_product=main_product, + main_product_chemical_IDs=main_product_chemical_IDs, + boiler=boiler, + add_EOL_GWP=True + ) + + # 打印结果 + print_microalgae_lca_results(lca) + + # 导出GWP分解数据 + export_gwp_breakdown_data(lca) + + # 导出完整结果到Excel + export_lca_results_to_excel(lca) + + # 生成详细报告 + generate_lca_report(lca) \ No newline at end of file diff --git a/biorefineries/microalgae/model_utils.py b/biorefineries/microalgae/model_utils.py new file mode 100644 index 00000000..bf7c1053 --- /dev/null +++ b/biorefineries/microalgae/model_utils.py @@ -0,0 +1,184 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Model utilities for microalgae biorefinery uncertainty analysis + +Based on succinic project's model_utils.py but adapted for microalgae system structure +""" +from pandas import DataFrame, read_excel +import chaospy as shape +import numpy as np +import biosteam as bst +from biosteam.evaluation import Model +from .system import microalgae_mcca_sys, microalgae_tea + +class MicroalgaeModel(bst.Model): + def __init__(self, system, metrics=None, specification=None, + parameters=None, retry_evaluation=True, exception_hook='warn', + namespace_dict={}): + Model.__init__(self, system=system, specification=specification, + parameters=parameters, retry_evaluation=retry_evaluation, exception_hook=exception_hook) + self.namespace_dict = namespace_dict + # Set metrics after initialization + if metrics is not None: + self.metrics = metrics + + def load_parameter_distributions(self, distributions, namespace_dict=None): + namespace_dict = namespace_dict or self.namespace_dict + + df = distributions + if type(df) is not DataFrame: + df = read_excel(distributions) + + create_function = self.create_function + param = self.parameter + + for i, row in df.iterrows(): + name = row['Parameter name'] + element = row['Element'] + kind = row['Kind'] + units = row['Units'] + baseline = row['Baseline'] + shape_data = row['Shape'] + lower, midpoint, upper = row['Lower'], row['Midpoint'], row['Upper'] + load_statements = str(row['Load statement']) + + D = None + if shape_data.lower() in ['triangular', 'triangle']: + D = shape.Triangle(lower, midpoint, upper) + elif shape_data.lower() in ['uniform']: + D = shape.Uniform(lower, upper) + + if D is not None: + param(name=name, + setter=create_function(load_statements, namespace_dict), + element=element, + kind=kind, + units=units, + baseline=baseline, + distribution=D) + + def create_function(self, code, namespace_dict): + def wrapper_fn(statement): + def f(x): + namespace_dict['x'] = x + exec(statement, namespace_dict) + return f + return wrapper_fn(code) + +def create_unit_groups(): + units_dict = {unit.ID: unit for unit in microalgae_mcca_sys.units} + unit_groups = [] + + # Area 1: Microalgae cultivation and harvesting + cultivation_units = [units_dict[uid] for uid in ['U101'] if uid in units_dict] + if cultivation_units: + unit_groups.append(bst.UnitGroup('Cultivation and harvesting', units=cultivation_units)) + + # Area 2: Pretreatment and hydrolysis + pretreatment_unit_ids = ['T201', 'P201', 'M201', 'P202', 'H201', 'R201', + 'T202', 'P203', 'R202', 'P204', 'H202', 'T203', + 'P205', 'T204', 'P206', 'T205', 'P207', 'S201', + 'M202', 'R203', 'H203', 'M203', 'R204', 'S202', 'P208'] + pretreatment_units = [units_dict[uid] for uid in pretreatment_unit_ids if uid in units_dict] + if pretreatment_units: + unit_groups.append(bst.UnitGroup('Pretreatment and hydrolysis', units=pretreatment_units)) + + # Area 3: Conversion + conversion_units = [units_dict[uid] for uid in ['H301', 'M301', 'T301', 'P301', 'R301', 'T302', 'S301'] if uid in units_dict] + if conversion_units: + unit_groups.append(bst.UnitGroup('Conversion', units=conversion_units)) + + # Area 4: Separation + separation_units = [units_dict[uid] for uid in ['M401', 'S402', 'D401', 'D402', 'D403', 'D404', 'D405'] if uid in units_dict] + if separation_units: + unit_groups.append(bst.UnitGroup('Separation', units=separation_units)) + + # Area 5: Waste treatment - Anaerobic digestion + waste_units = [units_dict[uid] for uid in ['M501', 'R501', 'M502', 'M503'] if uid in units_dict] + if waste_units: + unit_groups.append(bst.UnitGroup('Waste treatment and biogas', units=waste_units)) + + # Area 6: Storage + storage_unit_ids = ['T601', 'P601', 'T602', 'P602', 'T603', 'P603', 'T604', 'P604', 'T605', 'P605'] + storage_units = [units_dict[uid] for uid in storage_unit_ids if uid in units_dict] + if storage_units: + unit_groups.append(bst.UnitGroup('Storage', units=storage_units)) + + # Wastewater + wastewater_units = [units_dict[uid] for uid in ['WastewaterT'] if uid in units_dict] + remaining_units = [u for u in microalgae_mcca_sys.units + if u not in [unit for group in unit_groups for unit in group.units] + and u.ID not in ['WastewaterT']] + wastewater_units.extend(remaining_units) + if wastewater_units: + unit_groups.append(bst.UnitGroup('Wastewater treatment', units=wastewater_units)) + + # Boiler & turbogenerator + bt_units = [units_dict[uid] for uid in ['BT601'] if uid in units_dict] + if bt_units: + unit_groups.append(bst.UnitGroup('BT', units=bt_units)) + + # Heat Exchanger Network + hxn_units = [units_dict[uid] for uid in ['HXN601'] if uid in units_dict] + if hxn_units: + hxn_group = bst.UnitGroup('Heat exchange network', units=hxn_units) + hxn_group.filter_savings = False + unit_groups.append(hxn_group) + + # Other facilities + facility_units = [units_dict[uid] for uid in ['CT', 'PWC', 'ADP', 'CWP'] if uid in units_dict] + remaining_units = [u for u in microalgae_mcca_sys.units + if u not in [unit for group in unit_groups for unit in group.units] + and u.ID not in ['CT', 'PWC', 'ADP', 'CWP']] + facility_units.extend(remaining_units) + if facility_units: + unit_groups.append(bst.UnitGroup('Other facilities', units=facility_units)) + + # Fixed Operating Costs + foc_group = bst.UnitGroup('Fixed operating costs') + unit_groups.append(foc_group) + for ug in unit_groups: + ug.autofill_metrics(shorthand=False, + electricity_production=False, + electricity_consumption=True, + material_cost=True) + + # Special metric configurations + for ug in unit_groups: + if ug.name in ['Storage', 'Other facilities'] and ug.metrics: + for metric in ug.metrics: + if metric.name.lower() == 'material cost': + metric.getter = lambda: 0.0 + break + + # BT unit group + for ug in unit_groups: + if ug.name == 'BT' and ug.metrics: + for i, metric in enumerate(ug.metrics): + if metric.name.lower() == 'material cost': + ug.metrics[i] = bst.evaluation.Metric( + 'Material cost', + getter=lambda: microalgae_tea.utility_cost / microalgae_tea.operating_days / 24, + units='USD/hr', + element=None + ) + break + + # Fixed Operating Costs + for ug in unit_groups: + if ug.name == 'Fixed operating costs' and ug.metrics: + for i, metric in enumerate(ug.metrics): + if metric.name.lower() == 'material cost': + ug.metrics[i] = bst.evaluation.Metric( + 'Material cost', + getter=lambda: microalgae_tea.FOC / microalgae_tea.operating_days / 24, + units='USD/hr', + element=None + ) + break + + return unit_groups + +def get_unit_groups(): + return create_unit_groups() \ No newline at end of file diff --git a/biorefineries/microalgae/parameter_distributions.xlsx b/biorefineries/microalgae/parameter_distributions.xlsx new file mode 100644 index 00000000..05505cdc Binary files /dev/null and b/biorefineries/microalgae/parameter_distributions.xlsx differ diff --git a/biorefineries/microalgae/run_analysis.py b/biorefineries/microalgae/run_analysis.py new file mode 100644 index 00000000..12dfc8b7 --- /dev/null +++ b/biorefineries/microalgae/run_analysis.py @@ -0,0 +1,185 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Microalgae biorefinery analysis runner + +This script runs uncertainty and TRY analysis for the microalgae biorefinery. + +@author: Xingdong Shi +@version: 0.0.1 +""" + +import sys +import os +from datetime import datetime + +def run_uncertainty_analysis(): + """Run uncertainty analysis.""" + print("=" * 60) + print("MICROALGAE BIOREFINERY UNCERTAINTY ANALYSIS") + print("=" * 60) + + try: + from . import uncertainties + + print("Starting uncertainty analysis...") + results_dict, results = uncertainties.run_uncertainty_analysis() + + print("\nSaving results...") + uncertainties.save_results(results_dict, results) + + print("Uncertainty analysis completed successfully!") + return True + + except Exception as e: + print(f"Error in uncertainty analysis: {str(e)}") + import traceback + traceback.print_exc() + return False + +def run_try_analysis(): + """Run TRY analysis.""" + print("\n" + "=" * 60) + print("MICROALGAE BIOREFINERY TRY ANALYSIS") + print("=" * 60) + + try: + from . import TRY_analysis + + print("Starting TRY analysis...") + results_df = TRY_analysis.run_try_analysis() + + if results_df is not None: + print("\nAnalyzing results...") + analysis = TRY_analysis.analyze_try_results(results_df) + + print("Saving results...") + TRY_analysis.save_try_results(results_df, analysis) + + print("Running single parameter analysis...") + single_param_df = TRY_analysis.run_single_parameter_analysis() + + if single_param_df is not None: + timestamp = datetime.now().strftime('%Y%m%d_%H%M%S') + single_param_filename = f'microalgae_single_param_{timestamp}.xlsx' + single_param_df.to_excel(single_param_filename, index=False) + print(f'Single parameter results saved to {single_param_filename}') + + print("TRY analysis completed successfully!") + return True + else: + print("TRY analysis failed to generate results.") + return False + + except Exception as e: + print(f"Error in TRY analysis: {str(e)}") + import traceback + traceback.print_exc() + return False + +def run_quick_baseline(): + """Run a quick baseline analysis to verify system functionality.""" + print("\n" + "=" * 60) + print("MICROALGAE BIOREFINERY BASELINE CHECK") + print("=" * 60) + + try: + from . import system as microalgae_system + from . import lca + + print("Loading system...") + microalgae_sys = microalgae_system.microalgae_mcca_sys + microalgae_tea = microalgae_system.microalgae_tea + + print("Creating LCA object...") + from . import analysis_utils + microalgae_lca = analysis_utils.create_microalgae_lca_simple(microalgae_sys, microalgae_tea) + + print("Running baseline simulation...") + microalgae_sys.simulate() + + print("Calculating metrics...") + from . import analysis_utils + main_product = analysis_utils.get_main_product_stream(microalgae_sys) + mpsp = microalgae_tea.solve_price(main_product) if main_product else float('nan') + gwp = microalgae_lca.GWP + fec = microalgae_lca.FEC + + print(f"\nBaseline Results:") + print(f"MPSP: {mpsp:.3f} USD/kg") + print(f"GWP100a: {gwp:.3f} kg CO2-eq/kg") + print(f"FEC: {fec:.3f} MJ/kg") + + # Save baseline results + timestamp = datetime.now().strftime('%Y%m%d_%H%M%S') + baseline_filename = f'microalgae_baseline_{timestamp}.txt' + + with open(baseline_filename, 'w') as f: + f.write("Microalgae Biorefinery Baseline Results\n") + f.write("=" * 40 + "\n\n") + f.write(f"Timestamp: {datetime.now().strftime('%Y-%m-%d %H:%M:%S')}\n\n") + f.write(f"MPSP: {mpsp:.6f} USD/kg\n") + f.write(f"GWP100a: {gwp:.6f} kg CO2-eq/kg\n") + f.write(f"FEC: {fec:.6f} MJ/kg\n\n") + + # System information + f.write("System Information:\n") + f.write(f"Number of units: {len(microalgae_sys.units)}\n") + f.write(f"Number of streams: {len(microalgae_sys.streams)}\n") + f.write(f"Number of feeds: {len(microalgae_sys.feeds)}\n") + f.write(f"Number of products: {len(microalgae_sys.products)}\n") + + print(f"Baseline results saved to {baseline_filename}") + print("Baseline check completed successfully!") + return True + + except Exception as e: + print(f"Error in baseline check: {str(e)}") + import traceback + traceback.print_exc() + return False + +def main(): + """Main function to run all analyses.""" + print("MICROALGAE BIOREFINERY COMPREHENSIVE ANALYSIS") + print("=" * 60) + print(f"Start time: {datetime.now().strftime('%Y-%m-%d %H:%M:%S')}") + + # Run baseline check first + baseline_success = run_quick_baseline() + + if not baseline_success: + print("\nBaseline check failed. Stopping analysis.") + return + + # Ask user which analysis to run + print("\nWhich analysis would you like to run?") + print("1. Uncertainty analysis only") + print("2. TRY analysis only") + print("3. Both analyses") + print("4. Exit") + + choice = input("\nEnter your choice (1-4): ").strip() + + if choice == '1': + run_uncertainty_analysis() + elif choice == '2': + run_try_analysis() + elif choice == '3': + uncertainty_success = run_uncertainty_analysis() + if uncertainty_success: + run_try_analysis() + elif choice == '4': + print("Exiting...") + return + else: + print("Invalid choice. Running both analyses by default...") + uncertainty_success = run_uncertainty_analysis() + if uncertainty_success: + run_try_analysis() + + print(f"\nEnd time: {datetime.now().strftime('%Y-%m-%d %H:%M:%S')}") + print("Analysis completed!") + +if __name__ == '__main__': + main() \ No newline at end of file diff --git a/biorefineries/microalgae/sesitivity.py b/biorefineries/microalgae/sesitivity.py new file mode 100644 index 00000000..f9b1cf1b --- /dev/null +++ b/biorefineries/microalgae/sesitivity.py @@ -0,0 +1,776 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +BioSTEAM Sensitivity Analysis for Microalgae MCCA Biorefinery +""" + +import numpy as np +import pandas as pd +import matplotlib.pyplot as plt +import biosteam as bst +from biosteam.evaluation import Model, Metric +from chaospy import distributions as shape +import os +from matplotlib.colors import LinearSegmentedColormap +import seaborn as sns +from scipy.stats import spearmanr +import warnings +from microalgae.system import compute_labor_cost +from microalgae._chemicals import chems +from biosteam import main_flowsheet +from microalgae import system +from microalgae.system import create_microalgae_MCCA_production_sys +from microalgae.tea import create_tea + +# Suppress hot utility load warnings +warnings.filterwarnings("ignore", message="Hot utility load is negative.*", category=RuntimeWarning) + +# Set up plotting configuration +def setup_plotting_config(): + """Configure plotting environment""" + plt.rcParams['font.family'] = 'sans-serif' + plt.rcParams['font.sans-serif'] = ['Arial', 'DejaVu Sans'] + plt.rcParams['axes.unicode_minus'] = False + +# Initialize system +def initialize_system(): + """Initialize the biorefinery system and TEA""" + # Set up thermodynamics + bst.settings.set_thermo(chems) + main_flowsheet.clear() + + # Create system + from microalgae.system import create_microalgae_MCCA_production_sys + system = create_microalgae_MCCA_production_sys() + system.simulate() + + # Create TEA object + from microalgae.tea import create_tea + u = system.flowsheet.unit + dry_tpd = u.U101.ins[0].F_mass * 24 / 1000 # kg/h -> t/d + tea = create_tea( + system=system, + IRR=0.10, + duration=(2024, 2045), + depreciation='MACRS7', + income_tax=0.21, + operating_days=330, + lang_factor=None, + construction_schedule=(0.08, 0.60, 0.32), + startup_months=3, + startup_FOCfrac=1, + startup_salesfrac=0.5, + startup_VOCfrac=0.75, + WC_over_FCI=0.05, + finance_interest=0.08, + finance_years=10, + finance_fraction=0.4, + OSBL_units=(u.CT, u.CWP, u.ADP, u.PWC), + warehouse=0.04, + site_development=0.09, + additional_piping=0.045, + proratable_costs=0.10, + field_expenses=0.10, + construction=0.20, + contingency=0.10, + other_indirect_costs=0.10, + labor_cost=max(0.5e6, compute_labor_cost(dry_tpd)), + labor_burden=0.90, + property_insurance=0.007, + maintenance=0.03, + boiler_turbogenerator=None, + steam_power_depreciation='MACRS20' + ) + + return system, tea + +# Create evaluation model +def create_evaluation_model(microalgae_sys, microalgae_tea): + from biosteam.evaluation import Model, Metric + import chaospy as cp + + u = microalgae_sys.flowsheet.unit # 修正：使用flowsheet.unit而不是unit + s = microalgae_sys.flowsheet.stream # 修正：使用flowsheet.stream而不是stream + tea = microalgae_tea + + # 首先定义所有setter函数 + def set_microalgae_price(price): + u.U101.ins[0].price = price + + def set_caproic_acid_price(price): + s.caproic_acid_product.price = price + + def set_caproic_acid_yield_factor(factor): + u.R301.caproic_acid_yield_factor = factor + + def set_operating_days(days): + tea.operating_days = days + + def set_maintenance_factor(factor): + tea.maintenance = factor + + # 创建模型 + model = Model(microalgae_sys) + + # 定义指标 + model.metrics = [ + Metric('ROI', lambda: tea.ROI * 100, '%'), + Metric('MSP', lambda: tea.solve_price(s.caproic_acid_product), '$/kg'), + Metric('TCI', lambda: tea.TCI, '$'), + Metric('VOC', lambda: tea.VOC, '$'), + Metric('FOC', lambda: tea.FOC, '$') + ] + + # 之后再注册参数 - 修改为使用命名参数 + model.parameter(setter=set_microalgae_price, element=u.U101.ins[0], + name='Microalgae Price', kind='isolated', + units='$/kg', distribution=shape.Triangle(-0.2, -0.1, 0.1)) + + model.parameter(setter=set_caproic_acid_price, element=s.caproic_acid_product, + name='Caproic Acid Price', kind='isolated', + units='$/kg', distribution=shape.Triangle(2.0, 2.89, 4.0)) + + model.parameter(setter=set_caproic_acid_yield_factor, element=u.R301, + name='Caproic Acid Yield Factor', kind='isolated', + distribution=shape.Triangle(0.7, 1.0, 1.3)) + + model.parameter(setter=set_operating_days, element=tea, + name='Operating Days', kind='isolated', + distribution=shape.Triangle(300, 330, 350)) + + model.parameter(setter=set_maintenance_factor, element=tea, + name='Maintenance Factor', kind='isolated', + distribution=shape.Triangle(0.02, 0.03, 0.05)) + + # Print registered parameters to debug + print("Registered parameters:") + for i, param in enumerate(model.parameters): + print(f"{i+1}. {param.name} (index: {param.index})") + + return model + +# Univariate sensitivity analysis +def run_univariate_analysis(model, parameter_name, values=None, metrics=None): + """Run univariate sensitivity analysis""" + # Find parameter + parameter = None + for param in model.parameters: + if param.name == parameter_name: + parameter = param + break + + if parameter is None: + raise ValueError(f"Parameter '{parameter_name}' not found") + + # Generate value range + if values is None: + dist = parameter.distribution + if hasattr(dist, 'lower') and hasattr(dist, 'upper'): + values = np.linspace(dist.lower, dist.upper, 11) + else: + raise ValueError(f"Parameter '{parameter_name}' has no defined distribution range") + + print(f"Running univariate sensitivity analysis: {parameter_name}") + print(f"Parameter range: {min(values)} - {max(values)}") + + # Create sample points + baseline = model.get_baseline_sample() + results = [] + + # Evaluate each value + for value in values: + sample = baseline.copy() + sample[parameter.index] = value + # Directly call model for results + result = model(sample) + result[parameter.name] = value + results.append(result) + + # Combine results + df = pd.DataFrame(results) + + return df + +# Bivariate sensitivity analysis +def run_bivariate_analysis(model, param_name1, param_name2, metric_name='ROI', + values1=None, values2=None, save=True): + """Run bivariate sensitivity analysis""" + # Find matching parameters + param1 = param2 = None + for param in model.parameters: + if param.name == param_name1: + param1 = param + elif param.name == param_name2: + param2 = param + if param1 and param2: + break + + if not param1: + raise ValueError(f"Parameter '{param_name1}' not found") + if not param2: + raise ValueError(f"Parameter '{param_name2}' not found") + + # Set default value ranges + if values1 is None: + if hasattr(param1, 'distribution'): + dist = param1.distribution + if hasattr(dist, 'lower') and hasattr(dist, 'upper'): + values1 = np.linspace(dist.lower, dist.upper, 6) + + if values2 is None: + if hasattr(param2, 'distribution'): + dist = param2.distribution + if hasattr(dist, 'lower') and hasattr(dist, 'upper'): + values2 = np.linspace(dist.lower, dist.upper, 5) + + # Ensure values1 and values2 are 1D arrays + values1 = np.ravel(values1) + values2 = np.ravel(values2) + + print(f"Running bivariate sensitivity analysis: {param_name1} vs {param_name2}") + print(f"{param_name1} range: {min(values1)} - {max(values1)}") + print(f"{param_name2} range: {min(values2)} - {max(values2)}") + + # Get baseline sample + baseline = model.get_baseline_sample() + + # Create result grid + grid_data = np.zeros((len(values1), len(values2))) + + # Find correct metric index + metric_index = None + for i, metric in enumerate(model.metrics): + if metric.name == metric_name: + metric_index = i + break + + if metric_index is None: + raise ValueError(f"Metric '{metric_name}' not found. Available metrics: {[m.name for m in model.metrics]}") + + # Evaluate each parameter combination + for i, val1 in enumerate(values1): + for j, val2 in enumerate(values2): + sample = baseline.copy() + sample[param1.index] = val1 + sample[param2.index] = val2 + + # Directly evaluate model to get series result + result_series = model(sample) + # Use iloc instead of direct indexing + grid_data[i, j] = result_series.iloc[metric_index] if metric_index is not None else result_series[metric_name] + + # Create DataFrame for saving and plotting + result_df = pd.DataFrame(grid_data, index=values1, columns=values2) + + # Save results + if save: + os.makedirs('results', exist_ok=True) + safe_name1 = 'param1_' + str(hash(param_name1) % 10000) + safe_name2 = 'param2_' + str(hash(param_name2) % 10000) + safe_metric = 'metric_' + str(hash(metric_name) % 10000) + result_df.to_csv(f'results/bivariate_{safe_name1}_{safe_name2}_{safe_metric}.csv') + + # Create colormap + if metric_name == 'ROI': + cmap = 'RdYlGn' + center = 0 + else: + cmap = 'viridis' + center = None + + # Plot heatmap + plt.figure(figsize=(10, 8)) + ax = sns.heatmap(result_df, annot=True, cmap=cmap, center=center, + fmt=".2f", # Two decimal places + annot_kws={"size": 9}, # Adjust annotation size + cbar_kws={'label': f'{metric_name}'}) + + ax.set_title(f'Two-Variable Sensitivity Analysis: {metric_name}') + ax.set_xlabel(param_name2) + ax.set_ylabel(param_name1) + plt.tight_layout() + + # Save figure + if save: + os.makedirs('figures', exist_ok=True) + plt.savefig(f'figures/heatmap_{safe_name1}_{safe_name2}_{safe_metric}.png', dpi=300, bbox_inches='tight') + + return result_df + +# Monte Carlo analysis +def run_monte_carlo(model, n_samples=1000, sample_type='LHS', save=True): + """Run Monte Carlo uncertainty analysis using manual sample generation""" + print(f"Running Monte Carlo simulation ({n_samples} samples)...") + + # Get parameter distributions + param_distributions = {} + for param in model.parameters: + if hasattr(param, 'distribution'): + param_distributions[param.index] = param.distribution + + # Get metric name mapping + metric_names = {} + for i, metric in enumerate(model.metrics): + metric_names[i] = metric.name # Store position and metric name mapping + + # Get baseline sample as template + baseline = model.get_baseline_sample() + results = [] + + # Generate samples and evaluate each one + for i in range(n_samples): + if i % 100 == 0 and i > 0: + print(f"Completed {i}/{n_samples} samples") + + # Copy baseline sample + sample = baseline.copy() + + # Generate random values for each parameter + for param_name, distribution in param_distributions.items(): + if sample_type == 'LHS': + # Simplified Latin Hypercube Sampling + u = (i + 0.5) / n_samples # Ensure uniform coverage + if hasattr(distribution, 'ppf'): + sample[param_name] = distribution.ppf(u) + else: + # Fall back to regular random sampling + sample[param_name] = distribution.sample() + else: + # Regular random sampling + sample[param_name] = distribution.sample() + + # Evaluate sample + try: + # Get metric values + result_series = model(sample) + + # Create dictionary with parameter and metric values + result_dict = {} + + # Add parameter values + for param_name in param_distributions.keys(): + result_dict[param_name] = sample[param_name] + + # Add metric values using our defined names + for i, name in metric_names.items(): + result_dict[name] = result_series.iloc[i] + + results.append(result_dict) + except Exception as e: + print(f"Sample {i} evaluation failed: {str(e)}") + + # Combine results + df_results = pd.DataFrame(results) + + # Save results + if save: + os.makedirs('results', exist_ok=True) + df_results.to_csv('results/monte_carlo_results.csv', index=False) + + # Plot histograms - use correct column names + for metric_name in metric_names.values(): + if metric_name in df_results.columns: # Confirm column exists + plt.figure(figsize=(10, 6)) + + # Plot histogram and kernel density estimate + sns.histplot(df_results[metric_name], kde=True) + + # Add vertical lines for mean and median + mean_val = df_results[metric_name].mean() + median_val = df_results[metric_name].median() + plt.axvline(mean_val, color='r', linestyle='-', + label=f'Mean: {mean_val:.3f}') + plt.axvline(median_val, color='g', linestyle='--', + label=f'Median: {median_val:.3f}') + + plt.xlabel(metric_name) + plt.ylabel('Frequency') + plt.title(f'Monte Carlo Simulation: {metric_name} Distribution') + plt.legend() + + if save: + os.makedirs('figures', exist_ok=True) + plt.savefig(f'figures/monte_carlo_{metric_name}.png', dpi=300) + + return df_results + +# Run scenario analysis +def run_scenario_analysis(model, scenarios, save=True): + """Run predefined scenario analysis""" + print("Running scenario analysis...") + + # Get metric name mapping + metric_names = {} + for i, metric in enumerate(model.metrics): + metric_names[i] = metric.name + + # Create result containers + metrics = list(metric_names.values()) + results = {metric: [] for metric in metrics} + results['Scenario'] = [] + + # Get baseline sample as template + baseline = model.get_baseline_sample() + + for scenario_name, params in scenarios.items(): + results['Scenario'].append(scenario_name) + print(f"\nApplying scenario: {scenario_name}") + + # Display parameters + for param_name, value in params.items(): + print(f" {param_name} = {value}") + + # Create scenario sample + sample = baseline.copy() + + # Set scenario parameter values + for param in model.parameters: + if param.name in params: + sample[param.index] = params[param.name] + + # Directly evaluate model + try: + result_series = model(sample) + + # Add metric results + for i, metric_name in metric_names.items(): + value = result_series.iloc[i] + results[metric_name].append(value) + print(f" {metric_name}: {value}") + except Exception as e: + print(f"Scenario evaluation failed: {str(e)}") + # Add missing values + for metric_name in metric_names.values(): + results[metric_name].append(None) + + # Create result DataFrame + df_results = pd.DataFrame(results) + + # Save results + if save: + os.makedirs('results', exist_ok=True) + df_results.to_csv('results/scenario_analysis.csv', index=False) + + # Plot bar charts + for metric_name in metric_names.values(): + plt.figure(figsize=(10, 6)) + plt.bar(df_results['Scenario'], df_results[metric_name]) + plt.title(f'{metric_name} Across Scenarios') + plt.ylabel(f'{metric_name}') + plt.xticks(rotation=45) + plt.tight_layout() + + if save: + os.makedirs('figures', exist_ok=True) + plt.savefig(f'figures/scenario_{metric_name}.png', dpi=300, bbox_inches='tight') + + return df_results + +def plot_contour_from_heatmap(result_df, param_name1, param_name2, metric_name, + cmap='viridis', levels=15, save=True): + """Convert heatmap data to contour plot, similar to cane project style""" + fig, ax = plt.subplots(figsize=(10, 8)) + + # Get coordinate grid + x = result_df.columns.astype(float) + y = result_df.index.astype(float) + X, Y = np.meshgrid(x, y) + Z = result_df.values + + # Create custom colormap + if metric_name == 'ROI': + # ROI uses red-to-green gradient + cmap = LinearSegmentedColormap.from_list("ROI", ["#FF5555", "#FFFFFF", "#55FF55"]) + # Find zero point as center + vmin, vmax = Z.min(), Z.max() + if vmin < 0 < vmax: + vabs = max(abs(vmin), abs(vmax)) + norm = plt.Normalize(-vabs, vabs) + else: + norm = plt.Normalize(vmin, vmax) + else: + norm = None + + # Plot contour + cs = ax.contourf(X, Y, Z, levels=levels, cmap=cmap, norm=norm) + cs2 = ax.contour(X, Y, Z, levels=levels, colors='k', linewidths=0.5, alpha=0.3) + + # Add labels + ax.clabel(cs2, inline=True, fontsize=8, fmt='%.2f') + + # Add colorbar + cbar = plt.colorbar(cs, ax=ax) + cbar.set_label(metric_name) + + # Set labels + ax.set_xlabel(param_name2) + ax.set_ylabel(param_name1) + ax.set_title(f'Contour Map: {metric_name} as function of {param_name1} and {param_name2}') + + # Save figure + if save: + os.makedirs('figures', exist_ok=True) + safe_name1 = 'param1_' + str(hash(param_name1) % 10000) + safe_name2 = 'param2_' + str(hash(param_name2) % 10000) + safe_metric = 'metric_' + str(hash(metric_name) % 10000) + plt.savefig(f'figures/contour_{safe_name1}_{safe_name2}_{safe_metric}.png', + dpi=300, bbox_inches='tight') + + return fig, ax + +def plot_improved_sensitivity(df, parameter_name, metrics=None, save=True): + """Plot improved sensitivity curves similar to cane project style""" + if metrics is None: + metrics = [col for col in df.columns if col != parameter_name] + + # Create separate chart for each metric + for metric in metrics: + fig, ax = plt.subplots(figsize=(10, 6)) + + # Get data + x = df[parameter_name] + y = df[metric] + + # Plot curve with markers + ax.plot(x, y, '-o', lw=2.5, markersize=8, color='#3366CC', + markerfacecolor='white', markeredgecolor='#3366CC', markeredgewidth=2) + + # Add grid + ax.grid(True, linestyle='--', alpha=0.6) + + # Set title and labels + ax.set_title(f'Sensitivity: {metric} vs {parameter_name}', fontsize=14) + ax.set_xlabel(parameter_name, fontsize=12) + ax.set_ylabel(metric, fontsize=12) + + # Beautify axes + ax.spines['top'].set_visible(False) + ax.spines['right'].set_visible(False) + + # Save chart + if save: + os.makedirs('figures', exist_ok=True) + safe_param = 'param_' + str(hash(parameter_name) % 10000) + safe_metric = 'metric_' + str(hash(metric) % 10000) + plt.savefig(f'figures/sensitivity_{safe_param}_{safe_metric}_improved.png', + dpi=300, bbox_inches='tight') + + return fig, ax + +def plot_improved_monte_carlo(df, metric_name, save=True): + """Plot enhanced Monte Carlo distribution plot similar to cane project style""" + fig, ax = plt.subplots(figsize=(12, 7)) + + # Separate histogram and KDE curve plotting + # First plot histogram without KDE + sns.histplot(df[metric_name], kde=False, bins=25, ax=ax, + color='#3366CC', edgecolor='white', alpha=0.7) + + # Then plot KDE curve separately + sns.kdeplot(df[metric_name], ax=ax, color='red', lw=2) + + # Calculate and add statistics + mean_val = df[metric_name].mean() + median_val = df[metric_name].median() + std_val = df[metric_name].std() + + # Add vertical lines + ax.axvline(mean_val, color='red', linestyle='-', lw=2, + label=f'Mean: {mean_val:.3f}') + ax.axvline(median_val, color='green', linestyle='--', lw=2, + label=f'Median: {median_val:.3f}') + + # Add distribution area + p10 = df[metric_name].quantile(0.1) + p90 = df[metric_name].quantile(0.9) + ax.axvspan(p10, p90, alpha=0.2, color='grey', + label=f'10-90% CI: [{p10:.3f}, {p90:.3f}]') + + # Set title and labels + ax.set_title(f'Monte Carlo Distribution: {metric_name}', fontsize=14) + ax.set_xlabel(metric_name, fontsize=12) + ax.set_ylabel('Frequency', fontsize=12) + + # Add legend + ax.legend() + + # Save chart + if save: + os.makedirs('figures', exist_ok=True) + safe_metric = 'metric_' + str(hash(metric_name) % 10000) + plt.savefig(f'figures/monte_carlo_{safe_metric}_improved.png', + dpi=300, bbox_inches='tight') + + return fig, ax + +def plot_tornado(df, metric_name, parameters=None, n_parameters=10, save=True): + """Plot tornado diagram showing parameter influence on specific metric""" + # Calculate Spearman correlation coefficients + if parameters is None: + # Exclude metric columns, keep only parameter columns + parameters = [col for col in df.columns if col != metric_name] + + # Calculate correlation for each parameter with the metric + correlations = [] + for param in parameters: + corr, _ = spearmanr(df[param], df[metric_name]) + correlations.append((param, corr)) + + # Sort by correlation coefficient absolute value + correlations.sort(key=lambda x: abs(x[1]), reverse=True) + + # Keep only top n_parameters + if n_parameters < len(correlations): + correlations = correlations[:n_parameters] + + # Extract data for plotting + param_names, corr_values = zip(*correlations) + + # Plot tornado diagram + fig, ax = plt.subplots(figsize=(12, 8)) + + # Determine colors + colors = ['#3366CC' if c > 0 else '#CC3366' for c in corr_values] + + # Plot horizontal bars + bars = ax.barh(range(len(param_names)), corr_values, color=colors, + height=0.6, edgecolor='white', linewidth=1) + + # Add parameter name labels + ax.set_yticks(range(len(param_names))) + ax.set_yticklabels(param_names) + + # Set axes + ax.set_xlabel('Spearman Correlation Coefficient', fontsize=12) + ax.set_title(f'Parameter Importance for {metric_name}', fontsize=14) + + # Add vertical line at zero + ax.axvline(0, color='black', linestyle='-', linewidth=0.5) + + # Add grid + ax.grid(axis='x', linestyle='--', alpha=0.3) + + # Beautify + ax.spines['top'].set_visible(False) + ax.spines['right'].set_visible(False) + + # Add value labels for each bar + for i, bar in enumerate(bars): + width = bar.get_width() + label_x_pos = width + 0.01 if width >= 0 else width - 0.08 + ax.text(label_x_pos, bar.get_y() + bar.get_height()/2, + f'{width:.3f}', va='center') + + # Save chart + if save: + os.makedirs('figures', exist_ok=True) + safe_metric = 'metric_' + str(hash(metric_name) % 10000) + plt.savefig(f'figures/tornado_{safe_metric}.png', dpi=300, bbox_inches='tight') + + return fig, ax + +# Main program +if __name__ == "__main__": + print("Running BioSTEAM Sensitivity Analysis...") + + # Initialize system and model + bst.settings.set_thermo(chems) + main_flowsheet.clear() + microalgae_sys = create_microalgae_MCCA_production_sys() + u = microalgae_sys.flowsheet.unit + s = microalgae_sys.flowsheet.stream + + # Create TEA object + dry_tpd = u.U101.ins[0].F_mass * 24 / 1000 # kg/h -> t/d + microalgae_tea = create_tea( + system=microalgae_sys, + IRR=0.10, + duration=(2024, 2045), + depreciation='MACRS7', + income_tax=0.21, + operating_days=330, + lang_factor=None, + construction_schedule=(0.08, 0.60, 0.32), + startup_months=3, + startup_FOCfrac=1, + startup_salesfrac=0.5, + startup_VOCfrac=0.75, + WC_over_FCI=0.05, + finance_interest=0.08, + finance_years=10, + finance_fraction=0.4, + OSBL_units=(u.CT, u.CWP, u.ADP, u.PWC), + warehouse=0.04, + site_development=0.09, + additional_piping=0.045, + proratable_costs=0.10, + field_expenses=0.10, + construction=0.20, + contingency=0.10, + other_indirect_costs=0.10, + labor_cost=max(0.5e6, compute_labor_cost(dry_tpd)), + labor_burden=0.90, + property_insurance=0.007, + maintenance=0.03, + boiler_turbogenerator=None, + steam_power_depreciation='MACRS20' + ) + + # Create evaluation model + model = create_evaluation_model(microalgae_sys, microalgae_tea) + + # Define scenarios + scenarios = { + 'Optimistic': { + 'Microalgae price': -0.2, # Higher tipping fee income + 'Caproic acid price': 3.5, # Higher product price + 'Caproic acid yield factor': 1.2, # Higher than baseline yield + 'Operating days': 350, # More operating days + 'Maintenance factor': 0.02, # Lower maintenance cost + }, + 'Baseline': { + 'Microalgae price': -0.1, # Medium tipping fee income + 'Caproic acid price': 2.89, + 'Caproic acid yield factor': 1.0, + 'Operating days': 330, + 'Maintenance factor': 0.03, + }, + 'Pessimistic': { + 'Microalgae price': 0.0, # No tipping fee income + 'Caproic acid price': 2.0, # Lower product price + 'Caproic acid yield factor': 0.8, # Lower than baseline yield + 'Operating days': 300, # Fewer operating days + 'Maintenance factor': 0.05, # Higher maintenance cost + }, + 'Purchase': { + 'Microalgae price': 0.1, # Need to purchase feedstock + 'Caproic acid price': 3.2, # Higher product price + 'Caproic acid yield factor': 1.15, # Higher yield + 'Operating days': 330, + 'Maintenance factor': 0.03, + } + } + + # Single variable sensitivity analysis + print("\n1. Analyzing the effect of microalgae price on ROI and MSP...") + microalgae_price_results = run_univariate_analysis(model, 'Microalgae price') + plot_improved_sensitivity(microalgae_price_results, 'Microalgae price') + + # Bivariate analysis + print("\n2. Analyzing the effect of microalgae price and caproic acid price on ROI...") + heatmap_results = run_bivariate_analysis(model, 'Microalgae price', 'Caproic acid price', 'ROI') + plot_contour_from_heatmap(heatmap_results, 'Microalgae price', 'Caproic acid price', 'ROI') + + # Monte Carlo analysis + print("\n3. Running Monte Carlo simulation...") + monte_carlo_results = run_monte_carlo(model, n_samples=10000) + plot_improved_monte_carlo(monte_carlo_results, 'ROI') + + # Tornado diagram + plot_tornado(monte_carlo_results, 'ROI') + + # Scenario analysis + print("\n4. Running scenario analysis...") + scenario_results = run_scenario_analysis(model, scenarios) + + print("\nSensitivity analysis complete! Results saved in 'results' and 'figures' folders") \ No newline at end of file diff --git a/biorefineries/microalgae/streams.py b/biorefineries/microalgae/streams.py new file mode 100644 index 00000000..277a2220 --- /dev/null +++ b/biorefineries/microalgae/streams.py @@ -0,0 +1,30 @@ +# -*- coding: utf-8 -*- +""" +Created on Sat June 28 10:00:00 2025 + +Microalgae biorefinery to produce medium chain fatty acids +by anaerobic fermentation without external electron donor addition- feed streams + +References +---------- +[1] BioSTEAM Documentation: + https://biosteam.readthedocs.io/en/latest/tutorial/Creating_a_System.html +[2] Cortes-Peña et al., BioSTEAM: A Fast and Flexible Platform for the Design, + Simulation, and Techno-Economic Analysis of Biorefineries under Uncertainty. + ACS Sustainable Chem. Eng. 2020, 8 (8), 3302–3310. +[3] 3-Hydroxypropionic acid biorefineries project: + https://github.com/BioSTEAMDevelopmentGroup/Bioindustrial-Park/tree/master/biorefineries/HP +[4] Succinic projest + https://github.com/BioSTEAMDevelopmentGroup/Bioindustrial-Park/tree/master/biorefineries/succinic + +@author: Xingdong Shi +@version: 0.0.1 +""" + +from biosteam import stream_kwargs + +microalgae_feed = stream_kwargs('microalgae', + Microalgae=5000, + units='kg/hr', + price = 0.15 # this is a baseline value for TEA (cost for cultivation, harvesting and trasporting) +) diff --git a/biorefineries/microalgae/system.py b/biorefineries/microalgae/system.py new file mode 100644 index 00000000..25e10ed2 --- /dev/null +++ b/biorefineries/microalgae/system.py @@ -0,0 +1,430 @@ +# -*- coding: utf-8 -*- +""" +Created on Sat July 05 12:50:00 2025 + +Microalgae biorefinery to produce medium chain fatty acids +by anaerobic fermentation without external electron donor addition- system + +References +---------- +[1] BioSTEAM Documentation: + https://biosteam.readthedocs.io/en/latest/tutorial/Creating_a_System.html +[2] Cortes-Peña et al., BioSTEAM: A Fast and Flexible Platform for the Design, + Simulation, and Techno-Economic Analysis of Biorefineries under Uncertainty. + ACS Sustainable Chem. Eng. 2020, 8 (8), 3302–3310. +[3] 3-Hydroxypropionic acid biorefineries project: + https://github.com/BioSTEAMDevelopmentGroup/Bioindustrial-Park/tree/master/biorefineries/HP +[4] Succinic projest + https://github.com/BioSTEAMDevelopmentGroup/Bioindustrial-Park/tree/master/biorefineries/succinic + +@author: Xingdong Shi +@version: 0.0.8 +""" + +import biosteam as bst +import numpy as np +import thermosteam as tmo +from biosteam import Stream, SystemFactory +from biosteam.units import Pump, StorageTank, HXutility, Mixer, Splitter, MultiStageMixerSettlers +from biosteam.facilities import AirDistributionPackage, ProcessWaterCenter, CoolingTower, ChilledWaterPackage, HeatExchangerNetwork +from biosteam import main_flowsheet +from .units import ( + FeedstockPreprocessing, AcidPretreatmentReactor, Saccharification, SolidLiquidSeparation, MCCAFermentation, + NeutralizationTank, AnaerobicDigestion +) +from .utils import price +from ._chemicals import chems +from .tea import create_tea +from .streams import microalgae_feed +import warnings +# Filter out specific warnings +warnings.filterwarnings("ignore", message="phase equilibrium solution results in negative flow rates") +warnings.filterwarnings("ignore", message=".*has no defined Dortmund groups.*") +warnings.filterwarnings("ignore", message=".*has been replaced in registry") +warnings.filterwarnings("ignore", category=bst.exceptions.CostWarning) +warnings.filterwarnings("ignore", message=".*moisture.*is smaller than the desired.*") +warnings.filterwarnings("ignore", message=".*moisture of influent.*is smaller than the desired.*") + +# ------------------------------------------------------------------ +# Utility: labor cost scaling with plant size +# ------------------------------------------------------------------ +def compute_labor_cost(dry_tpd: float, + base_tpd: float = 2205.0, + base_cost: float = 3212962.0, + exponent: float = 0.2, + floor_tpd: float = 100.0, + floor_cost: float = 0.5e6) -> float: + if dry_tpd < floor_tpd: + return floor_cost + return base_cost * (dry_tpd / base_tpd) ** exponent + +# Set up the main flowsheet and thermodynamic environment +bst.settings.set_thermo(chems) +main_flowsheet.clear() +flowsheet = bst.Flowsheet('MCCA') +bst.main_flowsheet.set_flowsheet(flowsheet) + +# System settings +bst.System.default_converge_method = 'wegstein' +bst.System.default_maxiter = 2000 +bst.System.default_molar_tolerance = 1e-1 +bst.System.default_relative_molar_tolerance = 1e-1 # supersedes absolute tolerance +bst.System.strict_convergence = False # True => throw exception if system does not converge; False => continue with unconverged system + +@SystemFactory( + ID='Microalgae_MCCA_production', + ins=[dict(microalgae_feed, thermo=chems)], + outs=[ + #dict(ID='butanol_product', thermo=chems), + dict(ID='butyric_acid_product', thermo=chems), + dict(ID='caproic_acid_product', thermo=chems), + dict(ID='heptanoic_acid_product', thermo=chems), + dict(ID='caprylic_acid_product', thermo=chems)] + ) +def create_microalgae_MCCA_production_sys(ins, outs): + # Set the thermodynamic package explicitly + tmo.settings.set_thermo(chems) + + # Main feed and product + microalgae_feed, = ins + ( #butanol_product, + butyric_acid_product, + caproic_acid_product, + heptanoic_acid_product, + caprylic_acid_product) = outs + + # Calculate all required stream properties based on feed + microalgae_mass = microalgae_feed.F_mass + microalgae_water_mass = microalgae_mass / 0.04 # 4% solid loading + microalgae_water = Stream('microalgae_water', Water=microalgae_water_mass, units='kg/hr') + # H2SO4 for microalgae biomass hydrolysis + acid_loading = 1.47 # g H2SO4 / g microalgae + acid_purity = 0.93 + water_mass = microalgae_mass * (1 - 0.04) / 0.04 + pure_H2SO4 = microalgae_mass * acid_loading + acid_solution_mass = pure_H2SO4 / acid_purity + water_mass_acid = acid_solution_mass * (1 - acid_purity) + SulfuricAcid = Stream('sulfuricacid', H2SO4=pure_H2SO4, Water=water_mass_acid, units='kg/hr', price=price['SulfuricAcid']) + + # Store reference for later specification + _sulfuric_acid_stream = SulfuricAcid + # Ammonium Hydroxide for neutralization + h2so4_mol = pure_H2SO4 * 1000 / 98 # mol mass + nh4oh_mol = h2so4_mol * 0.1 # preadjustment + nh4oh_mass = nh4oh_mol * 35 / 1000 # mol mass to mass + ammonium_hydroxide = Stream('ammonium_hydroxide', NH4OH=nh4oh_mass, units='kg/hr', price=price['AmmoniumHydroxide']) + # Enzyme dosages + glucoamylase_mass = float(microalgae_mass * 0.0011) # ref from cron project + alpha_amylase_mass = float(microalgae_mass * 0.0082) # ref from cron project + glucoamylase = Stream('glucoamylase', GlucoAmylase=glucoamylase_mass, units='kg/hr', price=price['GlucoAmylase']) + alpha_amylase = Stream('alpha_amylase', AlphaAmylase=alpha_amylase_mass, units='kg/hr', price=price['AlphaAmylase']) + + # Store references for later specification + _glucoamylase_stream = glucoamylase + _alpha_amylase_stream = alpha_amylase + # NaOH for pH adjustment + naoh1_mass = microalgae_mass * 0.02 # pH 4.5 + naoh2_mass = microalgae_mass * 0.02 # pH 5.5 + total_naoh_mass = naoh1_mass + naoh2_mass + naoh = Stream('naoh', NaOH=total_naoh_mass, units='kg/hr', price=price['NaOH']) + # Yeast addition in MCCA fermentation + yeast_mass = microalgae_mass * 5 / 75 + fresh_yeast = Stream('fresh_yeast', Yeast=yeast_mass, units='kg/hr', price=price['Yeast']) + yeast_recycle = Stream('yeast_recycle', Yeast=0, units='kg/hr') + yeast_feed = bst.Stream() + # OleylAlcohol for extraction + fresh_oleylalcohol = Stream('fresh_oleylalcohol', OleylAlcohol=200, units='kg/hr', price=price['OleylAlcohol']) + oleylalcohol_recycle = Stream('oleylalcohol_recycle', OleylAlcohol=50, units='kg/hr') + oleylalcohol_feed = bst.Stream() + # Assign prices to product streams + #butanol_product.price = price['Butanol'] + butyric_acid_product.price = price['ButyricAcid'] + caproic_acid_product.price = price['CaproicAcid'] + heptanoic_acid_product.price = price['HeptanoicAcid'] + caprylic_acid_product.price = price['CaprylicAcid'] + + # ===================== + # Area 1: Microalgae process + # ===================== + U101 = FeedstockPreprocessing('U101', microalgae_feed, thermo=chems) + + # ===================== + # Area 2: Hydrolysis and saccharification + # ===================== + T201 = StorageTank('T201', SulfuricAcid) + P201 = Pump('P201', T201-0, P=5e5, pump_type='Default') + M201 = Mixer('M201', [P201-0, microalgae_water, U101-0]) + P202 = Pump('P202', M201-0) + H201 = HXutility('H201', P202-0, T=121+273.15) + R201 = AcidPretreatmentReactor('R201', H201-0) + T202 = StorageTank('T202', ammonium_hydroxide) + P203 = Pump('P203', T202-0) + R202 = NeutralizationTank('R202', [R201-0, P203-0]) + P204 = Pump('P204', R202-0) + H202 = HXutility('H202', P204-0, T=55+273.15) + T203 = StorageTank('T203', glucoamylase) + P205 = Pump('P205', T203-0, P=25e5, pump_type='Default', dP_design=24e5, ignore_NPSH=True) + T204 = StorageTank('T204', alpha_amylase) + P206 = Pump('P206', T204-0, P=25e5, pump_type='Default', dP_design=24e5, ignore_NPSH=True) + T205 = StorageTank('T205', naoh) + P207 = Pump('P207', T205-0, P=25e5, pump_type='Default', dP_design=24e5, ignore_NPSH=True) + S201 = Splitter('S201', P207-0, split=naoh1_mass/total_naoh_mass) + M202 = Mixer('M202', [H202-0, P205-0, S201-0]) + R203 = Saccharification('R203', M202-0) + H203 = HXutility('H203', R203-0, T=90+273.15) + M203 = Mixer('M203', [H203-0, P206-0, S201-1]) + R204 = Saccharification('R204', M203-0) + S202 = SolidLiquidSeparation('S202', R204-0) + P208 = Pump('P208', S202-0) + + # ===================== + # Area 3: Fermentation for MCCA production + # ===================== + H301 = HXutility('H301', P208-0, T=37+273.15) + + # Yeast mixing for fresh and recycle + M301 = Mixer('M301', [fresh_yeast, yeast_recycle], yeast_feed) + @M301.add_specification(run=True) + def adjust_fresh_yeast(): + total_yeast_needed = yeast_mass + recycle = yeast_recycle.imass['Yeast'] + fresh = max(total_yeast_needed - recycle, total_yeast_needed * 0.05) # 5% fresh yeast supply + fresh_yeast.imass['Yeast'] = fresh + + T301 = StorageTank('T301', yeast_feed) + P301 = Pump('P301', T301-0) + R301 = MCCAFermentation('R301', [H301-0, P301-0], ['', '', '', yeast_recycle], microalgae_mass_flow=microalgae_mass, titer = 2.003) + + # Add C6 yield factor specification + @R301.add_specification(run=True) + def set_C6_yield_factor(factor=1.0): + R301.caproic_acid_yield_factor = factor + + # Add C6 titer specification + @R301.add_specification(run=True) + def set_C6_titer(titer=2.003): + R301.titer = titer + + # Bind specification functions as methods for parameter loading + R301.set_C6_yield_factor = lambda x: setattr(R301, 'caproic_acid_yield_factor', x) + R301.set_C6_titer = lambda x: setattr(R301, 'titer', x) + + T302 = Mixer('T302', [R301-1]) + S301 = SolidLiquidSeparation('S301', R301-0) + + # ===================== + # Area 4: Product extraction + # ===================== + M401 = Mixer('M401', [fresh_oleylalcohol, oleylalcohol_recycle], oleylalcohol_feed) + @M401.add_specification(run=True) + def adjust_fresh_oleylalcohol(): + total_oleylalcohol = 200 + recycle = oleylalcohol_recycle.imass['OleylAlcohol'] + fresh = max(total_oleylalcohol - recycle, 1e-3) + fresh_oleylalcohol.imass['OleylAlcohol'] = fresh + IDs = ['Water', 'AceticAcid', 'PropionicAcid', 'ButyricAcid', 'ValericAcid', 'CaproicAcid','CaprylicAcid', 'HeptanoicAcid', 'Butanol', 'OleylAlcohol'] + K = np.array([1/5000, 0.24, 1.29, 5000/1, 13.58, 5000/1, 5000/1, 5000/1, 5000/1, 100000/1]) + S402 = MultiStageMixerSettlers( + 'S402', + partition_data={'K': K, 'IDs': IDs}, + N_stages=5, + ins=[S301-0, M401-0] + ) + + # Add extraction efficiency specification + original_K = K.copy() + @S402.add_specification(run=True) + def set_extraction_efficiency(efficiency=1.0): + S402.partition_data['K'] = original_K * efficiency + + # Bind as method for parameter loading + S402.set_extraction_efficiency = lambda x: setattr(S402, 'partition_data', {'K': original_K * x, 'IDs': IDs}) + + #D401 = ButanolDistillation('D401', S402-0) + #D401.check_LHK = False + D402 = bst.BinaryDistillation('D402', S402-0, LHK=('ButyricAcid', 'CaproicAcid'), + Lr=0.99, Hr=0.99, k=1.2, + partial_condenser=False, + is_divided=True) + #D402.check_LHK = False + D403 = bst.BinaryDistillation('D403', D402-1, LHK=('CaproicAcid', 'HeptanoicAcid'), + Lr=0.99, Hr=0.99, k=1.2, + partial_condenser=False, + is_divided=True + ) + #D403.check_LHK = False + D404 = bst.BinaryDistillation('D404', D403-1, LHK=('HeptanoicAcid', 'CaprylicAcid'), + Lr=0.99, Hr=0.99, k=1.2, + partial_condenser=False, + is_divided=True) + #D404.check_LHK = False + D405 = bst.BinaryDistillation('D405', D404-1, ['', oleylalcohol_recycle], LHK=('CaprylicAcid', 'OleylAlcohol'), + Lr=0.99, Hr=0.99, k=1.2, + partial_condenser=False, + is_divided=True, + product_specification_format='Recovery') + #D405.check_LHK = False + + # Add distillation efficiency specification + distillation_units = [D402, D403, D404, D405] + @D403.add_specification(run=True) # Use D403 as the representative unit + def set_distillation_efficiency(efficiency=1.0): + for unit in distillation_units: + unit.Lr = 0.99 * efficiency + unit.Hr = 0.99 * efficiency + + # Bind as method for parameter loading + D403.set_distillation_efficiency = lambda x: [setattr(unit, 'Lr', 0.99 * x) or setattr(unit, 'Hr', 0.99 * x) for unit in distillation_units] + + #Add acid loading factor specification to a suitable unit + @R201.add_specification(run=True) # Add to acid pretreatment reactor + def set_acid_loading_factor(factor=1.0): + pure_H2SO4_new = microalgae_mass * acid_loading * factor + acid_solution_mass_new = pure_H2SO4_new / acid_purity + water_mass_acid_new = acid_solution_mass_new * (1 - acid_purity) + _sulfuric_acid_stream.imass['H2SO4'] = pure_H2SO4_new + _sulfuric_acid_stream.imass['Water'] = water_mass_acid_new + + # Bind as method for parameter loading + def _set_acid_loading_factor(factor): + pure_H2SO4_new = microalgae_mass * acid_loading * factor + acid_solution_mass_new = pure_H2SO4_new / acid_purity + water_mass_acid_new = acid_solution_mass_new * (1 - acid_purity) + _sulfuric_acid_stream.imass['H2SO4'] = pure_H2SO4_new + _sulfuric_acid_stream.imass['Water'] = water_mass_acid_new + R201.set_acid_loading_factor = _set_acid_loading_factor + + # Add enzyme loading factor specification to a suitable unit + @R203.add_specification(run=True) # Add to first saccharification reactor + def set_enzyme_loading_factor(factor=1.0): + _glucoamylase_stream.imass['GlucoAmylase'] = microalgae_mass * 0.0011 * factor + _alpha_amylase_stream.imass['AlphaAmylase'] = microalgae_mass * 0.0082 * factor + + # Bind as method for parameter loading + def _set_enzyme_loading_factor(factor): + _glucoamylase_stream.imass['GlucoAmylase'] = microalgae_mass * 0.0011 * factor + _alpha_amylase_stream.imass['AlphaAmylase'] = microalgae_mass * 0.0082 * factor + R203.set_enzyme_loading_factor = _set_enzyme_loading_factor + + # ===================== + # Area 5: Waste reuse for biogas production + # ===================== + M501 = Mixer('M501', [S301-1, S202-1, S402-1]) + R501 = AnaerobicDigestion('R501', M501-0, microalgae_mass = microalgae_mass) + M503 = Mixer('M503', R501-1) + M502 = Mixer('M502', [R501-0, T302-0]) + + # ===================== + # Area 6: Facilities requirements + # ===================== + #T601 = StorageTank('T601', D401-0, tau=60.*24., V_wf=0.9, vessel_type='Floating roof', vessel_material='Stainless steel') + #P601 = Pump('P601', T601-0, butanol_product) + T602 = StorageTank('T602', D402-0, tau=30.*24., V_wf=0.9, vessel_type='Floating roof', vessel_material='Stainless steel') + P602 = Pump('P602', T602-0, butyric_acid_product) + T603 = StorageTank('T603', D403-0, tau=30.*24., V_wf=0.9, vessel_type='Floating roof', vessel_material='Stainless steel') + P603 = Pump('P603', T603-0, caproic_acid_product) + T604 = StorageTank('T604', D404-0, tau=60.*24., V_wf=0.9, vessel_type='Floating roof', vessel_material='Stainless steel') + P604 = Pump('P604', T604-0, heptanoic_acid_product) + T605 = StorageTank('T605', D405-0, tau=60.*24., V_wf=0.9, vessel_type='Floating roof', vessel_material='Stainless steel') + P605 = Pump('P605', T605-0, caprylic_acid_product) + CT = CoolingTower('CT') + HXN601 = HeatExchangerNetwork('HXN601', + #T_min_app=10, + #min_heat_util=2e6 + ) + PWC = ProcessWaterCenter('PWC') + ADP = AirDistributionPackage('ADP') + CWP = ChilledWaterPackage('CWP') + BT601 = bst.facilities.BoilerTurbogenerator('BT601', + ins=(R501-2, M502-0, '', '', '', ''), + satisfy_system_electricity_demand=True, + boiler_efficiency=0.9, + turbogenerator_efficiency=0.85) + + WastewaterT = bst.create_high_rate_wastewater_treatment_system('WastewaterT', + M503-0, # Use diluted wastewater stream + skip_IC=True, # Skip internal circulation to avoid division by zero + process_ID='6' # Use process ID 6 for unit numbering + ) + + +# ========================================== +# TEA Analysis +# ========================================== +# Create system and TEA objects at module level for import +u = flowsheet.unit +s = flowsheet.stream +microalgae_mcca_sys = create_microalgae_MCCA_production_sys() +microalgae_mcca_sys.simulate() + +# TEA analysis +# Dry biomass feed rate in ton per day (t/d) +dry_tpd = u.U101.ins[0].F_mass * 24 / 1000 # kg/h -> t/d +microalgae_tea = create_tea(system=microalgae_mcca_sys, IRR=0.10, duration=(2024, 2045), + depreciation='MACRS7', income_tax=0.21, + operating_days=330, + lang_factor= None, construction_schedule=(0.08, 0.60, 0.32), + startup_months=3, startup_FOCfrac=1, startup_salesfrac=0.5, + startup_VOCfrac=0.75, WC_over_FCI=0.05, + finance_interest=0.08, finance_years=10, finance_fraction=0.4, + OSBL_units=(u.CT, u.CWP, u.ADP, u.PWC, u.BT601), + warehouse=0.04, site_development=0.09, additional_piping=0.045, + proratable_costs=0.10, field_expenses=0.10, construction=0.20, + contingency=0.10, other_indirect_costs=0.10, + labor_cost=max(0.5e6, compute_labor_cost(dry_tpd)), + labor_burden=0.90, property_insurance=0.007, maintenance=0.03, boiler_turbogenerator=u.BT601, + steam_power_depreciation='MACRS20') + +if __name__ == '__main__': + microalgae_mcca_sys.diagram('cluster', format='png') + microalgae_mcca_sys.print() + print("\n===== Techno-Economic Analysis (TEA) Main Results =====") + # Use the system's main product stream directly for price calculation + caproic_acid_product = s.caproic_acid_product + if caproic_acid_product is not None: + if caproic_acid_product.price is None or caproic_acid_product.price == 0: + caproic_acid_product.price = 4.5 + price = microalgae_tea.solve_price(caproic_acid_product) + print(f"Caproic Acid Minimum Selling Price: {price:.2f} $/kg") + if caproic_acid_product.F_mass > 0 and caproic_acid_product.price > 0: + print("Caproic Acid Unit Production Cost:", microalgae_tea.production_costs([caproic_acid_product])) + print("NPV:", microalgae_tea.NPV) + print("TCI:", microalgae_tea.TCI) + print("FCI:", microalgae_tea.FCI) + print("DPI:", microalgae_tea.DPI) + print("TDC:", microalgae_tea.TDC) + print("FOC:", microalgae_tea.FOC) + print("VOC:", microalgae_tea.VOC) + print("AOC:", microalgae_tea.AOC) + print("ROI:", microalgae_tea.ROI) + print("PBP:", microalgae_tea.PBP) + print("Annual Depreciation:", microalgae_tea.annual_depreciation) + print("Sales:", microalgae_tea.sales) + print("Material Cost:", microalgae_tea.material_cost) + print("Utility Cost:", microalgae_tea.utility_cost) + print("CAPEX Table:\n", microalgae_tea.CAPEX_table()) + print("FOC Table:\n", microalgae_tea.FOC_table()) + print("Cashflow Table:\n", microalgae_tea.get_cashflow_table()) + + # Quick check: product flows in each units + # print("\n===== Stream Mass Flows for Each Unit (kg/hr) =====") + # for u_ in microalgae_mcca_sys.units: + # print(f"\n[{u_.ID} - {u_.__class__.__name__}]") + # for i, stream in enumerate(u_.ins): + # if stream: + # print(f" Inlet {i+1} ({stream.ID}):") + # for chem, flow in zip(stream.chemicals.IDs, stream.mass): + # if abs(flow) > 1e-6: + # print(f" {chem}: {flow:.2f} kg/hr") + # for i, stream in enumerate(u_.outs): + # if stream: + # print(f" Outlet {i+1} ({stream.ID}):") + # for chem, flow in zip(stream.chemicals.IDs, stream.mass): + # if abs(flow) > 1e-6: + # print(f" {chem}: {flow:.2f} kg/hr") + + # Quick check: product flows and prices + # for p in (#s.butanol_product, + # s.caproic_acid_product,s.heptanoic_acid_product, s.caprylic_acid_product, s.butyric_acid_product): + # print(f"{p.ID}: {p.F_mass:.2f} kg/h @ {p.price} $/kg") + diff --git a/biorefineries/microalgae/system_etanol.py b/biorefineries/microalgae/system_etanol.py new file mode 100644 index 00000000..fc291d79 --- /dev/null +++ b/biorefineries/microalgae/system_etanol.py @@ -0,0 +1,422 @@ +# -*- coding: utf-8 -*- +""" +Created on Sat July 05 12:50:00 2025 + +Microalgae biorefinery to produce medium chain fatty acids +by anaerobic fermentation without external electron donor addition- system + +References +---------- +[1] BioSTEAM Documentation: + https://biosteam.readthedocs.io/en/latest/tutorial/Creating_a_System.html +[2] Cortes-Peña et al., BioSTEAM: A Fast and Flexible Platform for the Design, + Simulation, and Techno-Economic Analysis of Biorefineries under Uncertainty. + ACS Sustainable Chem. Eng. 2020, 8 (8), 3302–3310. +[3] 3-Hydroxypropionic acid biorefineries project: + https://github.com/BioSTEAMDevelopmentGroup/Bioindustrial-Park/tree/master/biorefineries/HP +[4] Succinic projest + https://github.com/BioSTEAMDevelopmentGroup/Bioindustrial-Park/tree/master/biorefineries/succinic + +@author: Xingdong Shi +@version: 0.0.8 +""" + +import biosteam as bst +import numpy as np +import thermosteam as tmo +from biosteam import Stream, SystemFactory +from biosteam.units import Pump, StorageTank, HXutility, Mixer, Splitter, MultiStageMixerSettlers +from biosteam.facilities import AirDistributionPackage, ProcessWaterCenter, CoolingTower, ChilledWaterPackage, HeatExchangerNetwork +from biosteam import main_flowsheet +from .units import ( + FeedstockPreprocessing, AcidPretreatmentReactor, Saccharification, SolidLiquidSeparation, MCCAFermentation, + NeutralizationTank, AnaerobicDigestion +) +from .utils import price +from ._chemicals import chems +from .tea import create_tea +from .streams import microalgae_feed +import warnings +# Filter out specific warnings +# warnings.filterwarnings("ignore", message="phase equilibrium solution results in negative flow rates") +# warnings.filterwarnings("ignore", message=".*has no defined Dortmund groups.*") +# warnings.filterwarnings("ignore", message=".*has been replaced in registry") +# warnings.filterwarnings("ignore", category=bst.exceptions.CostWarning) +# warnings.filterwarnings("ignore", message=".*moisture.*is smaller than the desired.*") +# warnings.filterwarnings("ignore", message=".*moisture of influent.*is smaller than the desired.*") + +# ------------------------------------------------------------------ +# Utility: labor cost scaling with plant size +# ------------------------------------------------------------------ +def compute_labor_cost(dry_tpd: float, + base_tpd: float = 2205.0, + base_cost: float = 3212962.0, + exponent: float = 0.2, + floor_tpd: float = 100.0, + floor_cost: float = 0.5e6) -> float: + if dry_tpd < floor_tpd: + return floor_cost + return base_cost * (dry_tpd / base_tpd) ** exponent + +# Set up the main flowsheet and thermodynamic environment +bst.settings.set_thermo(chems) +main_flowsheet.clear() +flowsheet = bst.Flowsheet('MCCA_ethanol') +bst.main_flowsheet.set_flowsheet(flowsheet) + +# System settings +bst.System.default_converge_method = 'wegstein' +bst.System.default_maxiter = 2000 +bst.System.default_molar_tolerance = 1e-1 +bst.System.default_relative_molar_tolerance = 1e-1 # supersedes absolute tolerance +bst.System.strict_convergence = False # True => throw exception if system does not converge; False => continue with unconverged system + +@SystemFactory( + ID='Microalgae_MCCA_production_ethanol', + ins=[dict(microalgae_feed, thermo=chems)], + outs=[ + #dict(ID='butanol_product', thermo=chems), + dict(ID='butyric_acid_product', thermo=chems), + dict(ID='caproic_acid_product', thermo=chems), + dict(ID='heptanoic_acid_product', thermo=chems), + dict(ID='caprylic_acid_product', thermo=chems)] + ) +def create_microalgae_MCCA_production_sys_ethanol(ins, outs): + # Set the thermodynamic package explicitly + tmo.settings.set_thermo(chems) + + # Main feed and product + microalgae_feed, = ins + ( #butanol_product, + butyric_acid_product, + caproic_acid_product, + heptanoic_acid_product, + caprylic_acid_product) = outs + + # Calculate all required stream properties based on feed + microalgae_mass = microalgae_feed.F_mass + microalgae_water_mass = microalgae_mass / 0.04 # 4% solid loading + microalgae_water = Stream('microalgae_water', Water=microalgae_water_mass, units='kg/hr') + # H2SO4 for microalgae biomass hydrolysis + acid_loading = 1.47 # g H2SO4 / g microalgae + acid_purity = 0.93 + water_mass = microalgae_mass * (1 - 0.04) / 0.04 + pure_H2SO4 = microalgae_mass * acid_loading + acid_solution_mass = pure_H2SO4 / acid_purity + water_mass_acid = acid_solution_mass * (1 - acid_purity) + SulfuricAcid = Stream('sulfuricacid', H2SO4=pure_H2SO4, Water=water_mass_acid, units='kg/hr', price=price['SulfuricAcid']) + + # Store reference for later specification + _sulfuric_acid_stream = SulfuricAcid + # Ammonium Hydroxide for neutralization + h2so4_mol = pure_H2SO4 * 1000 / 98 # mol mass + nh4oh_mol = h2so4_mol * 0.1 # preadjustment + nh4oh_mass = nh4oh_mol * 35 / 1000 # mol mass to mass + ammonium_hydroxide = Stream('ammonium_hydroxide', NH4OH=nh4oh_mass, units='kg/hr', price=price['AmmoniumHydroxide']) + # Enzyme dosages + glucoamylase_mass = float(microalgae_mass * 0.0011) # ref from cron project + alpha_amylase_mass = float(microalgae_mass * 0.0082) # ref from cron project + glucoamylase = Stream('glucoamylase', GlucoAmylase=glucoamylase_mass, units='kg/hr', price=price['GlucoAmylase']) + alpha_amylase = Stream('alpha_amylase', AlphaAmylase=alpha_amylase_mass, units='kg/hr', price=price['AlphaAmylase']) + + # Store references for later specification + _glucoamylase_stream = glucoamylase + _alpha_amylase_stream = alpha_amylase + # NaOH for pH adjustment + naoh1_mass = microalgae_mass * 0.02 # pH 4.5 + naoh2_mass = microalgae_mass * 0.02 # pH 5.5 + total_naoh_mass = naoh1_mass + naoh2_mass + naoh = Stream('naoh', NaOH=total_naoh_mass, units='kg/hr', price=price['NaOH']) + # Ethanol addition in MCCA fermentation + # Accroding to the litterature, the dosgae of ethanol COD (2.09/g) would be ~2.0 times to microalgae COD (1.3475/g). + # Therefore, the etanol mass is 1.28 times of microalgae mass + # https://www.sciencedirect.com/science/article/pii/S0043135419309923?via%3Dihub + # https://www.sciencedirect.com/science/article/pii/S0045653521029465?via%3Dihub + + ethanol_mass = microalgae_mass * 1.28 + ethanol = Stream('ethanol', ethanol=ethanol_mass, units='kg/hr', price=price['Ethanol']) + # OleylAlcohol for extraction + fresh_oleylalcohol = Stream('fresh_oleylalcohol', OleylAlcohol=200, units='kg/hr', price=price['OleylAlcohol']) + oleylalcohol_recycle = Stream('oleylalcohol_recycle', OleylAlcohol=50, units='kg/hr') + oleylalcohol_feed = bst.Stream() + # Assign prices to product streams + #butanol_product.price = price['Butanol'] + butyric_acid_product.price = price['ButyricAcid'] + caproic_acid_product.price = price['CaproicAcid'] + heptanoic_acid_product.price = price['HeptanoicAcid'] + caprylic_acid_product.price = price['CaprylicAcid'] + + # ===================== + # Area 1: Microalgae process + # ===================== + U101 = FeedstockPreprocessing('U101', microalgae_feed, thermo=chems) + + # ===================== + # Area 2: Hydrolysis and saccharification + # ===================== + T201 = StorageTank('T201', SulfuricAcid) + P201 = Pump('P201', T201-0, P=5e5, pump_type='Default') + M201 = Mixer('M201', [P201-0, microalgae_water, U101-0]) + P202 = Pump('P202', M201-0) + H201 = HXutility('H201', P202-0, T=121+273.15) + R201 = AcidPretreatmentReactor('R201', H201-0) + T202 = StorageTank('T202', ammonium_hydroxide) + P203 = Pump('P203', T202-0) + R202 = NeutralizationTank('R202', [R201-0, P203-0]) + P204 = Pump('P204', R202-0) + H202 = HXutility('H202', P204-0, T=55+273.15) + T203 = StorageTank('T203', glucoamylase) + P205 = Pump('P205', T203-0, P=25e5, pump_type='Default', dP_design=24e5, ignore_NPSH=True) + T204 = StorageTank('T204', alpha_amylase) + P206 = Pump('P206', T204-0, P=25e5, pump_type='Default', dP_design=24e5, ignore_NPSH=True) + T205 = StorageTank('T205', naoh) + P207 = Pump('P207', T205-0, P=25e5, pump_type='Default', dP_design=24e5, ignore_NPSH=True) + S201 = Splitter('S201', P207-0, split=naoh1_mass/total_naoh_mass) + M202 = Mixer('M202', [H202-0, P205-0, S201-0]) + R203 = Saccharification('R203', M202-0) + H203 = HXutility('H203', R203-0, T=90+273.15) + M203 = Mixer('M203', [H203-0, P206-0, S201-1]) + R204 = Saccharification('R204', M203-0) + S202 = SolidLiquidSeparation('S202', R204-0) + P208 = Pump('P208', S202-0) + + # ===================== + # Area 3: Fermentation for MCCA production + # ===================== + H301 = HXutility('H301', P208-0, T=37+273.15) + T301 = StorageTank('T301', ethanol) + P301 = Pump('P301', T301-0) + R301 = MCCAFermentation('R301', [H301-0, P301-0], microalgae_mass_flow=microalgae_mass, titer = 2.003) + + # Add C6 yield factor specification + @R301.add_specification(run=True) + def set_C6_yield_factor(factor=1.0): + R301.caproic_acid_yield_factor = factor + + # Add C6 titer specification + @R301.add_specification(run=True) + def set_C6_titer(titer=2.003): + R301.titer = titer + + # Bind specification functions as methods for parameter loading + R301.set_C6_yield_factor = lambda x: setattr(R301, 'caproic_acid_yield_factor', x) + R301.set_C6_titer = lambda x: setattr(R301, 'titer', x) + + T302 = Mixer('T302', [R301-1]) + S301 = SolidLiquidSeparation('S301', R301-0) + + # ===================== + # Area 4: Product extraction + # ===================== + M401 = Mixer('M401', [fresh_oleylalcohol, oleylalcohol_recycle], oleylalcohol_feed) + @M401.add_specification(run=True) + def adjust_fresh_oleylalcohol(): + total_oleylalcohol = 200 + recycle = oleylalcohol_recycle.imass['OleylAlcohol'] + fresh = max(total_oleylalcohol - recycle, 1e-3) + fresh_oleylalcohol.imass['OleylAlcohol'] = fresh + IDs = ['Water', 'AceticAcid', 'PropionicAcid', 'ButyricAcid', 'ValericAcid', 'CaproicAcid','CaprylicAcid', 'HeptanoicAcid', 'Butanol', 'OleylAlcohol'] + K = np.array([1/5000, 0.24, 1.29, 5000/1, 13.58, 5000/1, 5000/1, 5000/1, 5000/1, 100000/1]) + S402 = MultiStageMixerSettlers( + 'S402', + partition_data={'K': K, 'IDs': IDs}, + N_stages=5, + ins=[S301-0, M401-0] + ) + + # Add extraction efficiency specification + original_K = K.copy() + @S402.add_specification(run=True) + def set_extraction_efficiency(efficiency=1.0): + S402.partition_data['K'] = original_K * efficiency + + # Bind as method for parameter loading + S402.set_extraction_efficiency = lambda x: setattr(S402, 'partition_data', {'K': original_K * x, 'IDs': IDs}) + + #D401 = ButanolDistillation('D401', S402-0) + #D401.check_LHK = False + D402 = bst.BinaryDistillation('D402', S402-0, LHK=('ButyricAcid', 'CaproicAcid'), + Lr=0.99, Hr=0.99, k=1.2, + partial_condenser=False, + is_divided=True) + #D402.check_LHK = False + D403 = bst.BinaryDistillation('D403', D402-1, LHK=('CaproicAcid', 'HeptanoicAcid'), + Lr=0.99, Hr=0.99, k=1.2, + partial_condenser=False, + is_divided=True + ) + #D403.check_LHK = False + D404 = bst.BinaryDistillation('D404', D403-1, LHK=('HeptanoicAcid', 'CaprylicAcid'), + Lr=0.99, Hr=0.99, k=1.2, + partial_condenser=False, + is_divided=True) + #D404.check_LHK = False + D405 = bst.BinaryDistillation('D405', D404-1, ['', oleylalcohol_recycle], LHK=('CaprylicAcid', 'OleylAlcohol'), + Lr=0.99, Hr=0.99, k=1.2, + partial_condenser=False, + is_divided=True, + product_specification_format='Recovery') + #D405.check_LHK = False + + # Add distillation efficiency specification + distillation_units = [D402, D403, D404, D405] + @D403.add_specification(run=True) # Use D403 as the representative unit + def set_distillation_efficiency(efficiency=1.0): + for unit in distillation_units: + unit.Lr = 0.99 * efficiency + unit.Hr = 0.99 * efficiency + + # Bind as method for parameter loading + D403.set_distillation_efficiency = lambda x: [setattr(unit, 'Lr', 0.99 * x) or setattr(unit, 'Hr', 0.99 * x) for unit in distillation_units] + + #Add acid loading factor specification to a suitable unit + @R201.add_specification(run=True) # Add to acid pretreatment reactor + def set_acid_loading_factor(factor=1.0): + pure_H2SO4_new = microalgae_mass * acid_loading * factor + acid_solution_mass_new = pure_H2SO4_new / acid_purity + water_mass_acid_new = acid_solution_mass_new * (1 - acid_purity) + _sulfuric_acid_stream.imass['H2SO4'] = pure_H2SO4_new + _sulfuric_acid_stream.imass['Water'] = water_mass_acid_new + + # Bind as method for parameter loading + def _set_acid_loading_factor(factor): + pure_H2SO4_new = microalgae_mass * acid_loading * factor + acid_solution_mass_new = pure_H2SO4_new / acid_purity + water_mass_acid_new = acid_solution_mass_new * (1 - acid_purity) + _sulfuric_acid_stream.imass['H2SO4'] = pure_H2SO4_new + _sulfuric_acid_stream.imass['Water'] = water_mass_acid_new + R201.set_acid_loading_factor = _set_acid_loading_factor + + # Add enzyme loading factor specification to a suitable unit + @R203.add_specification(run=True) # Add to first saccharification reactor + def set_enzyme_loading_factor(factor=1.0): + _glucoamylase_stream.imass['GlucoAmylase'] = microalgae_mass * 0.0011 * factor + _alpha_amylase_stream.imass['AlphaAmylase'] = microalgae_mass * 0.0082 * factor + + # Bind as method for parameter loading + def _set_enzyme_loading_factor(factor): + _glucoamylase_stream.imass['GlucoAmylase'] = microalgae_mass * 0.0011 * factor + _alpha_amylase_stream.imass['AlphaAmylase'] = microalgae_mass * 0.0082 * factor + R203.set_enzyme_loading_factor = _set_enzyme_loading_factor + + # ===================== + # Area 5: Waste reuse for biogas production + # ===================== + M501 = Mixer('M501', [S301-1, S202-1, S402-1]) + R501 = AnaerobicDigestion('R501', M501-0, microalgae_mass = microalgae_mass) + M503 = Mixer('M503', R501-1) + M502 = Mixer('M502', [R501-0, T302-0]) + + # ===================== + # Area 6: Facilities requirements + # ===================== + #T601 = StorageTank('T601', D401-0, tau=60.*24., V_wf=0.9, vessel_type='Floating roof', vessel_material='Stainless steel') + #P601 = Pump('P601', T601-0, butanol_product) + T602 = StorageTank('T602', D402-0, tau=30.*24., V_wf=0.9, vessel_type='Floating roof', vessel_material='Stainless steel') + P602 = Pump('P602', T602-0, butyric_acid_product) + T603 = StorageTank('T603', D403-0, tau=30.*24., V_wf=0.9, vessel_type='Floating roof', vessel_material='Stainless steel') + P603 = Pump('P603', T603-0, caproic_acid_product) + T604 = StorageTank('T604', D404-0, tau=60.*24., V_wf=0.9, vessel_type='Floating roof', vessel_material='Stainless steel') + P604 = Pump('P604', T604-0, heptanoic_acid_product) + T605 = StorageTank('T605', D405-0, tau=60.*24., V_wf=0.9, vessel_type='Floating roof', vessel_material='Stainless steel') + P605 = Pump('P605', T605-0, caprylic_acid_product) + CT = CoolingTower('CT') + HXN601 = HeatExchangerNetwork('HXN601', + #T_min_app=10, + #min_heat_util=2e6 + ) + PWC = ProcessWaterCenter('PWC') + ADP = AirDistributionPackage('ADP') + CWP = ChilledWaterPackage('CWP') + BT601 = bst.facilities.BoilerTurbogenerator('BT601', + ins=(R501-2, M502-0, '', '', '', ''), + satisfy_system_electricity_demand=True, + boiler_efficiency=0.9, + turbogenerator_efficiency=0.85) + + WastewaterT = bst.create_high_rate_wastewater_treatment_system('WastewaterT', + M503-0, # Use diluted wastewater stream + skip_IC=True, # Skip internal circulation to avoid division by zero + process_ID='6' # Use process ID 6 for unit numbering + ) + +# ========================================== +# TEA Analysis +# ========================================== +# Create system and TEA objects at module level for import +u = flowsheet.unit +s = flowsheet.stream +microalgae_mcca_sys_ethanol = create_microalgae_MCCA_production_sys_ethanol() +microalgae_mcca_sys_ethanol.simulate() + +# TEA analysis +# Dry biomass feed rate in ton per day (t/d) +dry_tpd = u.U101.ins[0].F_mass * 24 / 1000 # kg/h -> t/d +microalgae_tea_ethanol = create_tea(system=microalgae_mcca_sys_ethanol, IRR=0.10, duration=(2024, 2045), + depreciation='MACRS7', income_tax=0.21, + operating_days=330, + lang_factor= None, construction_schedule=(0.08, 0.60, 0.32), + startup_months=3, startup_FOCfrac=1, startup_salesfrac=0.5, + startup_VOCfrac=0.75, WC_over_FCI=0.05, + finance_interest=0.08, finance_years=10, finance_fraction=0.4, + OSBL_units=(u.CT, u.CWP, u.ADP, u.PWC, u.BT601), + warehouse=0.04, site_development=0.09, additional_piping=0.045, + proratable_costs=0.10, field_expenses=0.10, construction=0.20, + contingency=0.10, other_indirect_costs=0.10, + labor_cost=max(0.5e6, compute_labor_cost(dry_tpd)), + labor_burden=0.90, property_insurance=0.007, maintenance=0.03, boiler_turbogenerator=u.BT601, + steam_power_depreciation='MACRS20') + +if __name__ == '__main__': + microalgae_mcca_sys_ethanol.diagram('cluster', format='png') + microalgae_mcca_sys_ethanol.print() + print("\n===== Techno-Economic Analysis (TEA) Main Results =====") + # Use the system's main product stream directly for price calculation + caproic_acid_product = s.caproic_acid_product + if caproic_acid_product is not None: + if caproic_acid_product.price is None or caproic_acid_product.price == 0: + caproic_acid_product.price = 4.5 + price = microalgae_tea_ethanol.solve_price(caproic_acid_product) + print(f"Caproic Acid Minimum Selling Price: {price:.2f} $/kg") + if caproic_acid_product.F_mass > 0 and caproic_acid_product.price > 0: + print("Caproic Acid Unit Production Cost:", microalgae_tea_ethanol.production_costs([caproic_acid_product])) + print("NPV:", microalgae_tea_ethanol.NPV) + print("TCI:", microalgae_tea_ethanol.TCI) + print("FCI:", microalgae_tea_ethanol.FCI) + print("DPI:", microalgae_tea_ethanol.DPI) + print("TDC:", microalgae_tea_ethanol.TDC) + print("FOC:", microalgae_tea_ethanol.FOC) + print("VOC:", microalgae_tea_ethanol.VOC) + print("AOC:", microalgae_tea_ethanol.AOC) + print("ROI:", microalgae_tea_ethanol.ROI) + print("PBP:", microalgae_tea_ethanol.PBP) + print("Annual Depreciation:", microalgae_tea_ethanol.annual_depreciation) + print("Sales:", microalgae_tea_ethanol.sales) + print("Material Cost:", microalgae_tea_ethanol.material_cost) + print("Utility Cost:", microalgae_tea_ethanol.utility_cost) + print("CAPEX Table:\n", microalgae_tea_ethanol.CAPEX_table()) + print("FOC Table:\n", microalgae_tea_ethanol.FOC_table()) + print("Cashflow Table:\n", microalgae_tea_ethanol.get_cashflow_table()) + + # Quick check: product flows in each units + # print("\n===== Stream Mass Flows for Each Unit (kg/hr) =====") + # for u_ in microalgae_mcca_sys.units: + # print(f"\n[{u_.ID} - {u_.__class__.__name__}]") + # for i, stream in enumerate(u_.ins): + # if stream: + # print(f" Inlet {i+1} ({stream.ID}):") + # for chem, flow in zip(stream.chemicals.IDs, stream.mass): + # if abs(flow) > 1e-6: + # print(f" {chem}: {flow:.2f} kg/hr") + # for i, stream in enumerate(u_.outs): + # if stream: + # print(f" Outlet {i+1} ({stream.ID}):") + # for chem, flow in zip(stream.chemicals.IDs, stream.mass): + # if abs(flow) > 1e-6: + # print(f" {chem}: {flow:.2f} kg/hr") + + # Quick check: product flows and prices + # for p in (#s.butanol_product, + # s.caproic_acid_product,s.heptanoic_acid_product, s.caprylic_acid_product, s.butyric_acid_product): + # print(f"{p.ID}: {p.F_mass:.2f} kg/h @ {p.price} $/kg") + diff --git a/biorefineries/microalgae/system_no_yeast.py b/biorefineries/microalgae/system_no_yeast.py new file mode 100644 index 00000000..c289e496 --- /dev/null +++ b/biorefineries/microalgae/system_no_yeast.py @@ -0,0 +1,415 @@ +# -*- coding: utf-8 -*- +""" +Created on Sat July 05 12:50:00 2025 + +Microalgae biorefinery to produce medium chain fatty acids +by anaerobic fermentation without external electron donor addition- system + +References +---------- +[1] BioSTEAM Documentation: + https://biosteam.readthedocs.io/en/latest/tutorial/Creating_a_System.html +[2] Cortes-Peña et al., BioSTEAM: A Fast and Flexible Platform for the Design, + Simulation, and Techno-Economic Analysis of Biorefineries under Uncertainty. + ACS Sustainable Chem. Eng. 2020, 8 (8), 3302–3310. +[3] 3-Hydroxypropionic acid biorefineries project: + https://github.com/BioSTEAMDevelopmentGroup/Bioindustrial-Park/tree/master/biorefineries/HP +[4] Succinic projest + https://github.com/BioSTEAMDevelopmentGroup/Bioindustrial-Park/tree/master/biorefineries/succinic + +@author: Xingdong Shi +@version: 0.0.8 +""" + +import biosteam as bst +import numpy as np +import thermosteam as tmo +from biosteam import Stream, SystemFactory +from biosteam.units import Pump, StorageTank, HXutility, Mixer, Splitter, MultiStageMixerSettlers +from biosteam.facilities import AirDistributionPackage, ProcessWaterCenter, CoolingTower, ChilledWaterPackage, HeatExchangerNetwork +from biosteam import main_flowsheet +from .units import ( + FeedstockPreprocessing, AcidPretreatmentReactor, Saccharification, SolidLiquidSeparation, MCCAFermentation_no_yeast, + NeutralizationTank, AnaerobicDigestion +) +from .utils import price +from ._chemicals import chems +from .tea import create_tea +from .streams import microalgae_feed +import warnings +# Filter out specific warnings +# warnings.filterwarnings("ignore", message="phase equilibrium solution results in negative flow rates") +# warnings.filterwarnings("ignore", message=".*has no defined Dortmund groups.*") +# warnings.filterwarnings("ignore", message=".*has been replaced in registry") +# warnings.filterwarnings("ignore", category=bst.exceptions.CostWarning) +# warnings.filterwarnings("ignore", message=".*moisture.*is smaller than the desired.*") +# warnings.filterwarnings("ignore", message=".*moisture of influent.*is smaller than the desired.*") + +# ------------------------------------------------------------------ +# Utility: labor cost scaling with plant size +# ------------------------------------------------------------------ +def compute_labor_cost(dry_tpd: float, + base_tpd: float = 2205.0, + base_cost: float = 3212962.0, + exponent: float = 0.2, + floor_tpd: float = 100.0, + floor_cost: float = 0.5e6) -> float: + if dry_tpd < floor_tpd: + return floor_cost + return base_cost * (dry_tpd / base_tpd) ** exponent + +# Set up the main flowsheet and thermodynamic environment +bst.settings.set_thermo(chems) +main_flowsheet.clear() +flowsheet = bst.Flowsheet('MCCA') +bst.main_flowsheet.set_flowsheet(flowsheet) + +# System settings +bst.System.default_converge_method = 'wegstein' +bst.System.default_maxiter = 2000 +bst.System.default_molar_tolerance = 1e-1 +bst.System.default_relative_molar_tolerance = 1e-1 # supersedes absolute tolerance +bst.System.strict_convergence = False # True => throw exception if system does not converge; False => continue with unconverged system + +@SystemFactory( + ID='Microalgae_MCCA_production_no_yeast', + ins=[dict(microalgae_feed, thermo=chems)], + outs=[ + #dict(ID='butanol_product', thermo=chems), + dict(ID='butyric_acid_product', thermo=chems), + dict(ID='caproic_acid_product', thermo=chems), + dict(ID='heptanoic_acid_product', thermo=chems), + dict(ID='caprylic_acid_product', thermo=chems)] + ) +def create_microalgae_MCCA_production_no_yeast_sys(ins, outs): + # Set the thermodynamic package explicitly + tmo.settings.set_thermo(chems) + + # Main feed and product + microalgae_feed, = ins + ( #butanol_product, + butyric_acid_product, + caproic_acid_product, + heptanoic_acid_product, + caprylic_acid_product) = outs + + # Calculate all required stream properties based on feed + microalgae_mass = microalgae_feed.F_mass + microalgae_water_mass = microalgae_mass / 0.04 # 4% solid loading + microalgae_water = Stream('microalgae_water', Water=microalgae_water_mass, units='kg/hr') + # H2SO4 for microalgae biomass hydrolysis + acid_loading = 1.47 # g H2SO4 / g microalgae + acid_purity = 0.93 + water_mass = microalgae_mass * (1 - 0.04) / 0.04 + pure_H2SO4 = microalgae_mass * acid_loading + acid_solution_mass = pure_H2SO4 / acid_purity + water_mass_acid = acid_solution_mass * (1 - acid_purity) + SulfuricAcid = Stream('sulfuricacid', H2SO4=pure_H2SO4, Water=water_mass_acid, units='kg/hr', price=price['SulfuricAcid']) + + # Store reference for later specification + _sulfuric_acid_stream = SulfuricAcid + # Ammonium Hydroxide for neutralization + h2so4_mol = pure_H2SO4 * 1000 / 98 # mol mass + nh4oh_mol = h2so4_mol * 0.1 # preadjustment + nh4oh_mass = nh4oh_mol * 35 / 1000 # mol mass to mass + ammonium_hydroxide = Stream('ammonium_hydroxide', NH4OH=nh4oh_mass, units='kg/hr', price=price['AmmoniumHydroxide']) + # Enzyme dosages + glucoamylase_mass = float(microalgae_mass * 0.0011) # ref from cron project + alpha_amylase_mass = float(microalgae_mass * 0.0082) # ref from cron project + glucoamylase = Stream('glucoamylase', GlucoAmylase=glucoamylase_mass, units='kg/hr', price=price['GlucoAmylase']) + alpha_amylase = Stream('alpha_amylase', AlphaAmylase=alpha_amylase_mass, units='kg/hr', price=price['AlphaAmylase']) + + # Store references for later specification + _glucoamylase_stream = glucoamylase + _alpha_amylase_stream = alpha_amylase + # NaOH for pH adjustment + naoh1_mass = microalgae_mass * 0.02 # pH 4.5 + naoh2_mass = microalgae_mass * 0.02 # pH 5.5 + total_naoh_mass = naoh1_mass + naoh2_mass + naoh = Stream('naoh', NaOH=total_naoh_mass, units='kg/hr', price=price['NaOH']) + # OleylAlcohol for extraction + fresh_oleylalcohol = Stream('fresh_oleylalcohol', OleylAlcohol=200, units='kg/hr', price=price['OleylAlcohol']) + oleylalcohol_recycle = Stream('oleylalcohol_recycle', OleylAlcohol=50, units='kg/hr') + oleylalcohol_feed = bst.Stream() + # Assign prices to product streams + #butanol_product.price = price['Butanol'] + butyric_acid_product.price = price['ButyricAcid'] + caproic_acid_product.price = price['CaproicAcid'] + heptanoic_acid_product.price = price['HeptanoicAcid'] + caprylic_acid_product.price = price['CaprylicAcid'] + + # ===================== + # Area 1: Microalgae process + # ===================== + U101 = FeedstockPreprocessing('U101', microalgae_feed, thermo=chems) + + # ===================== + # Area 2: Hydrolysis and saccharification + # ===================== + T201 = StorageTank('T201', SulfuricAcid) + P201 = Pump('P201', T201-0, P=5e5, pump_type='Default') + M201 = Mixer('M201', [P201-0, microalgae_water, U101-0]) + P202 = Pump('P202', M201-0) + H201 = HXutility('H201', P202-0, T=121+273.15) + R201 = AcidPretreatmentReactor('R201', H201-0) + T202 = StorageTank('T202', ammonium_hydroxide) + P203 = Pump('P203', T202-0) + R202 = NeutralizationTank('R202', [R201-0, P203-0]) + P204 = Pump('P204', R202-0) + H202 = HXutility('H202', P204-0, T=55+273.15) + T203 = StorageTank('T203', glucoamylase) + P205 = Pump('P205', T203-0, P=25e5, pump_type='Default', dP_design=24e5, ignore_NPSH=True) + T204 = StorageTank('T204', alpha_amylase) + P206 = Pump('P206', T204-0, P=25e5, pump_type='Default', dP_design=24e5, ignore_NPSH=True) + T205 = StorageTank('T205', naoh) + P207 = Pump('P207', T205-0, P=25e5, pump_type='Default', dP_design=24e5, ignore_NPSH=True) + S201 = Splitter('S201', P207-0, split=naoh1_mass/total_naoh_mass) + M202 = Mixer('M202', [H202-0, P205-0, S201-0]) + R203 = Saccharification('R203', M202-0) + H203 = HXutility('H203', R203-0, T=90+273.15) + M203 = Mixer('M203', [H203-0, P206-0, S201-1]) + R204 = Saccharification('R204', M203-0) + S202 = SolidLiquidSeparation('S202', R204-0) + P208 = Pump('P208', S202-0) + + # ===================== + # Area 3: Fermentation for MCCA production + # ===================== + H301 = HXutility('H301', P208-0, T=37+273.15) + T301 = StorageTank('T301', H301-0) + P301 = Pump('P301', T301-0) + R301 = MCCAFermentation_no_yeast('R301', [H301-0, P301-0], microalgae_mass_flow=microalgae_mass, titer = 1.208) + + # Add C6 yield factor specification + @R301.add_specification(run=True) + def set_C6_yield_factor(factor=1.0): + R301.caproic_acid_yield_factor = factor + + # Add C6 titer specification + @R301.add_specification(run=True) + def set_C6_titer(titer=2.003): + R301.titer = titer + + # Bind specification functions as methods for parameter loading + R301.set_C6_yield_factor = lambda x: setattr(R301, 'caproic_acid_yield_factor', x) + R301.set_C6_titer = lambda x: setattr(R301, 'titer', x) + + T302 = Mixer('T302', [R301-1]) + S301 = SolidLiquidSeparation('S301', R301-0) + + # ===================== + # Area 4: Product extraction + # ===================== + M401 = Mixer('M401', [fresh_oleylalcohol, oleylalcohol_recycle], oleylalcohol_feed) + @M401.add_specification(run=True) + def adjust_fresh_oleylalcohol(): + total_oleylalcohol = 200 + recycle = oleylalcohol_recycle.imass['OleylAlcohol'] + fresh = max(total_oleylalcohol - recycle, 1e-3) + fresh_oleylalcohol.imass['OleylAlcohol'] = fresh + IDs = ['Water', 'AceticAcid', 'PropionicAcid', 'ButyricAcid', 'ValericAcid', 'CaproicAcid','CaprylicAcid', 'HeptanoicAcid', 'Butanol', 'OleylAlcohol'] + K = np.array([1/5000, 0.24, 1.29, 5000/1, 13.58, 5000/1, 5000/1, 5000/1, 5000/1, 100000/1]) + S402 = MultiStageMixerSettlers( + 'S402', + partition_data={'K': K, 'IDs': IDs}, + N_stages=5, + ins=[S301-0, M401-0] + ) + + # Add extraction efficiency specification + original_K = K.copy() + @S402.add_specification(run=True) + def set_extraction_efficiency(efficiency=1.0): + S402.partition_data['K'] = original_K * efficiency + + # Bind as method for parameter loading + S402.set_extraction_efficiency = lambda x: setattr(S402, 'partition_data', {'K': original_K * x, 'IDs': IDs}) + + #D401 = ButanolDistillation('D401', S402-0) + #D401.check_LHK = False + D402 = bst.BinaryDistillation('D402', S402-0, LHK=('ButyricAcid', 'CaproicAcid'), + Lr=0.99, Hr=0.99, k=1.2, + partial_condenser=False, + is_divided=True) + #D402.check_LHK = False + D403 = bst.BinaryDistillation('D403', D402-1, LHK=('CaproicAcid', 'HeptanoicAcid'), + Lr=0.99, Hr=0.99, k=1.2, + partial_condenser=False, + is_divided=True + ) + #D403.check_LHK = False + D404 = bst.BinaryDistillation('D404', D403-1, LHK=('HeptanoicAcid', 'CaprylicAcid'), + Lr=0.99, Hr=0.99, k=1.2, + partial_condenser=False, + is_divided=True) + #D404.check_LHK = False + D405 = bst.BinaryDistillation('D405', D404-1, ['', oleylalcohol_recycle], LHK=('CaprylicAcid', 'OleylAlcohol'), + Lr=0.99, Hr=0.99, k=1.2, + partial_condenser=False, + is_divided=True, + product_specification_format='Recovery') + #D405.check_LHK = False + + # Add distillation efficiency specification + distillation_units = [D402, D403, D404, D405] + @D403.add_specification(run=True) # Use D403 as the representative unit + def set_distillation_efficiency(efficiency=1.0): + for unit in distillation_units: + unit.Lr = 0.99 * efficiency + unit.Hr = 0.99 * efficiency + + # Bind as method for parameter loading + D403.set_distillation_efficiency = lambda x: [setattr(unit, 'Lr', 0.99 * x) or setattr(unit, 'Hr', 0.99 * x) for unit in distillation_units] + + #Add acid loading factor specification to a suitable unit + @R201.add_specification(run=True) # Add to acid pretreatment reactor + def set_acid_loading_factor(factor=1.0): + pure_H2SO4_new = microalgae_mass * acid_loading * factor + acid_solution_mass_new = pure_H2SO4_new / acid_purity + water_mass_acid_new = acid_solution_mass_new * (1 - acid_purity) + _sulfuric_acid_stream.imass['H2SO4'] = pure_H2SO4_new + _sulfuric_acid_stream.imass['Water'] = water_mass_acid_new + + # Bind as method for parameter loading + def _set_acid_loading_factor(factor): + pure_H2SO4_new = microalgae_mass * acid_loading * factor + acid_solution_mass_new = pure_H2SO4_new / acid_purity + water_mass_acid_new = acid_solution_mass_new * (1 - acid_purity) + _sulfuric_acid_stream.imass['H2SO4'] = pure_H2SO4_new + _sulfuric_acid_stream.imass['Water'] = water_mass_acid_new + R201.set_acid_loading_factor = _set_acid_loading_factor + + # Add enzyme loading factor specification to a suitable unit + @R203.add_specification(run=True) # Add to first saccharification reactor + def set_enzyme_loading_factor(factor=1.0): + _glucoamylase_stream.imass['GlucoAmylase'] = microalgae_mass * 0.0011 * factor + _alpha_amylase_stream.imass['AlphaAmylase'] = microalgae_mass * 0.0082 * factor + + # Bind as method for parameter loading + def _set_enzyme_loading_factor(factor): + _glucoamylase_stream.imass['GlucoAmylase'] = microalgae_mass * 0.0011 * factor + _alpha_amylase_stream.imass['AlphaAmylase'] = microalgae_mass * 0.0082 * factor + R203.set_enzyme_loading_factor = _set_enzyme_loading_factor + + # ===================== + # Area 5: Waste reuse for biogas production + # ===================== + M501 = Mixer('M501', [S301-1, S202-1, S402-1]) + R501 = AnaerobicDigestion('R501', M501-0, microalgae_mass = microalgae_mass) + M503 = Mixer('M503', R501-1) + M502 = Mixer('M502', [R501-0, T302-0]) + + # ===================== + # Area 6: Facilities requirements + # ===================== + #T601 = StorageTank('T601', D401-0, tau=60.*24., V_wf=0.9, vessel_type='Floating roof', vessel_material='Stainless steel') + #P601 = Pump('P601', T601-0, butanol_product) + T602 = StorageTank('T602', D402-0, tau=30.*24., V_wf=0.9, vessel_type='Floating roof', vessel_material='Stainless steel') + P602 = Pump('P602', T602-0, butyric_acid_product) + T603 = StorageTank('T603', D403-0, tau=30.*24., V_wf=0.9, vessel_type='Floating roof', vessel_material='Stainless steel') + P603 = Pump('P603', T603-0, caproic_acid_product) + T604 = StorageTank('T604', D404-0, tau=60.*24., V_wf=0.9, vessel_type='Floating roof', vessel_material='Stainless steel') + P604 = Pump('P604', T604-0, heptanoic_acid_product) + T605 = StorageTank('T605', D405-0, tau=60.*24., V_wf=0.9, vessel_type='Floating roof', vessel_material='Stainless steel') + P605 = Pump('P605', T605-0, caprylic_acid_product) + CT = CoolingTower('CT') + HXN601 = HeatExchangerNetwork('HXN601', + #T_min_app=10, + #min_heat_util=2e6 + ) + PWC = ProcessWaterCenter('PWC') + ADP = AirDistributionPackage('ADP') + CWP = ChilledWaterPackage('CWP') + BT601 = bst.facilities.BoilerTurbogenerator('BT601', + ins=(R501-2, M502-0, '', '', '', ''), + satisfy_system_electricity_demand=True, + boiler_efficiency=0.9, + turbogenerator_efficiency=0.85) + + WastewaterT = bst.create_high_rate_wastewater_treatment_system('WastewaterT', + M503-0, # Use diluted wastewater stream + skip_IC=True, # Skip internal circulation to avoid division by zero + process_ID='6' # Use process ID 6 for unit numbering + ) + + +# ========================================== +# TEA Analysis +# ========================================== +# Create system and TEA objects at module level for import +u = flowsheet.unit +s = flowsheet.stream +microalgae_mcca_sys_no_yeast = create_microalgae_MCCA_production_no_yeast_sys() +microalgae_mcca_sys_no_yeast.simulate() + +# TEA analysis +# Dry biomass feed rate in ton per day (t/d) +dry_tpd = u.U101.ins[0].F_mass * 24 / 1000 # kg/h -> t/d +microalgae_tea_no_yeast = create_tea(system=microalgae_mcca_sys_no_yeast, IRR=0.10, duration=(2024, 2045), + depreciation='MACRS7', income_tax=0.21, + operating_days=330, + lang_factor= None, construction_schedule=(0.08, 0.60, 0.32), + startup_months=3, startup_FOCfrac=1, startup_salesfrac=0.5, + startup_VOCfrac=0.75, WC_over_FCI=0.05, + finance_interest=0.08, finance_years=10, finance_fraction=0.4, + OSBL_units=(u.CT, u.CWP, u.ADP, u.PWC, u.BT601), + warehouse=0.04, site_development=0.09, additional_piping=0.045, + proratable_costs=0.10, field_expenses=0.10, construction=0.20, + contingency=0.10, other_indirect_costs=0.10, + labor_cost=max(0.5e6, compute_labor_cost(dry_tpd)), + labor_burden=0.90, property_insurance=0.007, maintenance=0.03, boiler_turbogenerator=u.BT601, + steam_power_depreciation='MACRS20') + +if __name__ == '__main__': + microalgae_mcca_sys_no_yeast.diagram('cluster', format='png') + microalgae_mcca_sys_no_yeast.print() + print("\n===== Techno-Economic Analysis (TEA) Main Results =====") + # Use the system's main product stream directly for price calculation + caproic_acid_product = s.caproic_acid_product + if caproic_acid_product is not None: + if caproic_acid_product.price is None or caproic_acid_product.price == 0: + caproic_acid_product.price = 4.5 + price = microalgae_tea_no_yeast.solve_price(caproic_acid_product) + print(f"Caproic Acid Minimum Selling Price: {price:.2f} $/kg") + if caproic_acid_product.F_mass > 0 and caproic_acid_product.price > 0: + print("Caproic Acid Unit Production Cost:", microalgae_tea_no_yeast.production_costs([caproic_acid_product])) + print("NPV:", microalgae_tea_no_yeast.NPV) + print("TCI:", microalgae_tea_no_yeast.TCI) + print("FCI:", microalgae_tea_no_yeast.FCI) + print("DPI:", microalgae_tea_no_yeast.DPI) + print("TDC:", microalgae_tea_no_yeast.TDC) + print("FOC:", microalgae_tea_no_yeast.FOC) + print("VOC:", microalgae_tea_no_yeast.VOC) + print("AOC:", microalgae_tea_no_yeast.AOC) + print("ROI:", microalgae_tea_no_yeast.ROI) + print("PBP:", microalgae_tea_no_yeast.PBP) + print("Annual Depreciation:", microalgae_tea_no_yeast.annual_depreciation) + print("Sales:", microalgae_tea_no_yeast.sales) + print("Material Cost:", microalgae_tea_no_yeast.material_cost) + print("Utility Cost:", microalgae_tea_no_yeast.utility_cost) + print("CAPEX Table:\n", microalgae_tea_no_yeast.CAPEX_table()) + print("FOC Table:\n", microalgae_tea_no_yeast.FOC_table()) + print("Cashflow Table:\n", microalgae_tea_no_yeast.get_cashflow_table()) + + # Quick check: product flows in each units + # print("\n===== Stream Mass Flows for Each Unit (kg/hr) =====") + # for u_ in microalgae_mcca_sys.units: + # print(f"\n[{u_.ID} - {u_.__class__.__name__}]") + # for i, stream in enumerate(u_.ins): + # if stream: + # print(f" Inlet {i+1} ({stream.ID}):") + # for chem, flow in zip(stream.chemicals.IDs, stream.mass): + # if abs(flow) > 1e-6: + # print(f" {chem}: {flow:.2f} kg/hr") + # for i, stream in enumerate(u_.outs): + # if stream: + # print(f" Outlet {i+1} ({stream.ID}):") + # for chem, flow in zip(stream.chemicals.IDs, stream.mass): + # if abs(flow) > 1e-6: + # print(f" {chem}: {flow:.2f} kg/hr") + + # Quick check: product flows and prices + # for p in (#s.butanol_product, + # s.caproic_acid_product,s.heptanoic_acid_product, s.caprylic_acid_product, s.butyric_acid_product): + # print(f"{p.ID}: {p.F_mass:.2f} kg/h @ {p.price} $/kg") + diff --git a/biorefineries/microalgae/tea.py b/biorefineries/microalgae/tea.py new file mode 100644 index 00000000..17f65dc3 --- /dev/null +++ b/biorefineries/microalgae/tea.py @@ -0,0 +1,27 @@ +# -*- coding: utf-8 -*- +""" +Created on Sat July 20 17:50:00 2025 + +Microalgae biorefinery to produce medium chain fatty acids +by anaerobic fermentation without external electron donor addition- TEA + +References +---------- +[1] BioSTEAM Documentation: + https://biosteam.readthedocs.io/en/latest/tutorial/Creating_a_System.html +[2] Cortes-Peña et al., BioSTEAM: A Fast and Flexible Platform for the Design, + Simulation, and Techno-Economic Analysis of Biorefineries under Uncertainty. + ACS Sustainable Chem. Eng. 2020, 8 (8), 3302–3310. +[3] 3-Hydroxypropionic acid biorefineries project: + https://github.com/BioSTEAMDevelopmentGroup/Bioindustrial-Park/tree/master/biorefineries/HP +[4] Succinic projest + https://github.com/BioSTEAMDevelopmentGroup/Bioindustrial-Park/tree/master/biorefineries/succinic + +@author: Xingdong Shi +@version: 0.0.8 +""" + +import biosteam as bst +from biorefineries.tea.cellulosic_ethanol_tea import CellulosicEthanolTEA as create_tea + + diff --git a/biorefineries/microalgae/tea_breakdown.py b/biorefineries/microalgae/tea_breakdown.py new file mode 100644 index 00000000..6667af12 --- /dev/null +++ b/biorefineries/microalgae/tea_breakdown.py @@ -0,0 +1,204 @@ +# -*- coding: utf-8 -*- +""" +Microalgae TEA Breakdown Analysis + +@author: Xingdong Shi +""" + +import pandas as pd +import numpy as np +import biosteam as bst +from .model_utils import get_unit_groups +from .system import microalgae_mcca_sys, microalgae_tea + +class BTElectricityRevenueConfig: + def __init__(self, + electricity_price=0.07, # USD/kWh + include_in_material_cost=True, + include_in_operating_cost=False): + self.electricity_price = electricity_price + self.include_in_material_cost = include_in_material_cost + self.include_in_operating_cost = include_in_operating_cost + +bt_revenue_config = BTElectricityRevenueConfig() + +def create_tea_breakdown_data(unit_groups=None, tea=microalgae_tea, print_output=False, fractions=False): + unit_groups = unit_groups or get_unit_groups() + system_electricity_revenue = calculate_system_electricity_revenue() + first_metrics_group = None + for ug in unit_groups: + if hasattr(ug, 'metrics') and ug.metrics: + first_metrics_group = ug + break + metric_breakdowns = {metric.name: {} for metric in first_metrics_group.metrics} + for ug in unit_groups: + if hasattr(ug, 'metrics') and ug.metrics: + for metric in ug.metrics: + denominator = 1.0 + + if fractions: + metric_name = metric.name + if 'cost' in metric_name.lower() or 'installed' in metric_name.lower(): + denominator = tea.installed_equipment_cost / 1e6 + elif 'cooling' in metric_name.lower(): + denominator = 0 + for unit in microalgae_mcca_sys.units: + if hasattr(unit, 'heat_utilities'): + for hu in unit.heat_utilities: + if hasattr(hu, 'duty') and hu.duty > 0: + denominator += hu.duty / 1e6 + elif 'heating' in metric_name.lower(): + denominator = 0 + for unit in microalgae_mcca_sys.units: + if hasattr(unit, 'heat_utilities'): + for hu in unit.heat_utilities: + if hasattr(hu, 'duty') and hu.duty < 0: + denominator -= hu.duty / 1e6 + elif 'electricity' in metric_name.lower() or 'power' in metric_name.lower(): + if 'consumption' in metric_name.lower(): + denominator = microalgae_mcca_sys.power_utility.consumption / 1e3 + elif 'production' in metric_name.lower(): + denominator = microalgae_mcca_sys.power_utility.production / 1e3 + elif 'material' in metric_name.lower(): + denominator = tea.material_cost / tea.operating_hours + metric_value = metric() if callable(metric) else 0 + if ug.name == 'BT': + should_add_revenue = ( + (bt_revenue_config.include_in_material_cost and 'material' in metric.name.lower()) or + (bt_revenue_config.include_in_operating_cost and ('operating' in metric.name.lower() or 'variable' in metric.name.lower())) + ) + if should_add_revenue: + if fractions: + electricity_revenue_to_add = system_electricity_revenue + else: + electricity_revenue_to_add = system_electricity_revenue / tea.operating_hours + metric_value += electricity_revenue_to_add + if ug.name == 'Storage': + continue + elif ug.name == 'Other facilities': + storage_value = 0 + for storage_ug in unit_groups: + if storage_ug.name == 'Storage': + if hasattr(storage_ug, 'metrics') and storage_ug.metrics: + storage_metric = storage_ug.metrics[ug.metrics.index(metric)] if ug.metrics.index(metric) < len(storage_ug.metrics) else None + storage_value = storage_metric() if callable(storage_metric) else 0 + break + + combined_name = 'Storage and ' + ug.name + metric_breakdowns[metric.name][combined_name] = (metric_value + storage_value) / denominator + else: + metric_breakdowns[metric.name][ug.name] = metric_value / denominator + if print_output and first_metrics_group: + for metric in first_metrics_group.metrics: + print(f"\n\n----- {metric.name} ({metric.units}) -----") + metric_breakdown = metric_breakdowns.get(metric.name, {}) + for group_name, value in metric_breakdown.items(): + print(f"{group_name}: {value:.3f}") + return metric_breakdowns + +def calculate_system_electricity_revenue(): + power_breakdown = get_system_power_breakdown() + net_power_demand = power_breakdown['net_electricity'] # kW + operating_hours = microalgae_tea.operating_hours + + if net_power_demand < 0: + surplus_electricity = abs(net_power_demand) # kW + annual_electricity_revenue = surplus_electricity * bt_revenue_config.electricity_price * operating_hours + return -annual_electricity_revenue + else: + electricity_cost = net_power_demand * bt_revenue_config.electricity_price * operating_hours + return electricity_cost + +def create_tea_breakdown_dataframe(unit_groups=None, fraction=True, + scale_fractions_to_positive_values=False): + unit_groups = unit_groups or get_unit_groups() + df = bst.UnitGroup.df_from_groups( + unit_groups, + fraction=fraction, + scale_fractions_to_positive_values=scale_fractions_to_positive_values + ) + + system_electricity_revenue = calculate_system_electricity_revenue() + if system_electricity_revenue != 0 and 'BT' in df.index: + if 'Material cost' in df.columns: + if fraction: + total_annual_material_cost = microalgae_tea.material_cost + revenue_fraction = system_electricity_revenue / total_annual_material_cost + df.loc['BT', 'Material cost'] += revenue_fraction + else: + revenue_hourly = system_electricity_revenue / microalgae_tea.operating_hours + df.loc['BT', 'Material cost'] += revenue_hourly + return df + +def get_cost_breakdown_by_category(tea): + breakdown = {} + breakdown['Installed equipment cost'] = tea.installed_equipment_cost / 1e6 + breakdown['ISBL installed equipment cost'] = getattr(tea, 'ISBL_installed_equipment_cost', 0) / 1e6 + breakdown['OSBL installed equipment cost'] = getattr(tea, 'OSBL_installed_equipment_cost', 0) / 1e6 + breakdown['Fixed operating cost'] = tea.FOC / 1e6 + breakdown['Variable operating cost'] = tea.VOC / 1e6 + breakdown['Material cost'] = tea.material_cost / 1e6 + breakdown['Utility cost'] = tea.utility_cost / 1e6 + breakdown['Total capital investment'] = tea.TCI / 1e6 + breakdown['Fixed capital investment'] = tea.FCI / 1e6 + breakdown['Total depreciable capital'] = tea.TDC / 1e6 + breakdown['Direct permanent investment'] = tea.DPI / 1e6 + return breakdown + +def get_system_power_breakdown(): + power_breakdown = { + 'total_consumption': 0, + 'total_production': 0, + 'net_electricity': 0, + 'unit_details': {} + } + + for unit in microalgae_mcca_sys.units: + if hasattr(unit, 'power_utility') and unit.power_utility: + pu = unit.power_utility + consumption = pu.consumption + production = pu.production + if consumption != 0 or production != 0: + power_breakdown['unit_details'][unit.ID] = { + 'consumption': consumption, + 'production': production, + 'net': consumption - production + } + power_breakdown['total_consumption'] += consumption + power_breakdown['total_production'] += production + power_breakdown['net_electricity'] = (power_breakdown['total_consumption'] - + power_breakdown['total_production']) + return power_breakdown + +def analyze_microalgae_tea_breakdown(): + microalgae_mcca_sys.simulate() + system_electricity_revenue = calculate_system_electricity_revenue() + power_breakdown = get_system_power_breakdown() + unit_groups = get_unit_groups() + breakdown_data = create_tea_breakdown_data(unit_groups, microalgae_tea, print_output=True) + df_breakdown = create_tea_breakdown_dataframe( + unit_groups, + fraction=True, + scale_fractions_to_positive_values=False + ) + df_breakdown = create_tea_breakdown_dataframe( + unit_groups, + fraction=True, + scale_fractions_to_positive_values=False + ) + cost_breakdown = get_cost_breakdown_by_category(microalgae_tea) + results = { + 'breakdown_data': breakdown_data, + 'breakdown_dataframe': df_breakdown, + 'cost_breakdown': cost_breakdown, + 'unit_groups': unit_groups, + 'tea_system': microalgae_mcca_sys, + 'power_breakdown': power_breakdown, + 'system_electricity_revenue': system_electricity_revenue + } + return results + +if __name__ == "__main__": + analyze_microalgae_tea_breakdown() + + diff --git a/biorefineries/microalgae/uncertainties.py b/biorefineries/microalgae/uncertainties.py new file mode 100644 index 00000000..d4ea2a92 --- /dev/null +++ b/biorefineries/microalgae/uncertainties.py @@ -0,0 +1,494 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +# Bioindustrial-Park: BioSTEAM's Premier Biorefinery Models and Results +# Copyright (C) 2021-, Sarang Bhagwat +# +# This module is under the UIUC open-source license. See +# github.com/BioSTEAMDevelopmentGroup/biosteam/blob/master/LICENSE.txt +# for license details. + +Uncertainty analysis for microalgae biorefinery + +Based on succinic project but adapted for microalgae system structure +""" + +from warnings import filterwarnings +filterwarnings('ignore') +import numpy as np +import pandas as pd +import biosteam as bst +# import contourplots +from .system import microalgae_mcca_sys, microalgae_tea +from .lca import create_microalgae_lca +from .model_utils import MicroalgaeModel +from biosteam.evaluation import Metric +from datetime import datetime +from biosteam.utils import TicToc +import os + +chdir = os.chdir +microalgae_filepath = os.path.dirname(__file__) +microalgae_results_filepath = os.path.join(microalgae_filepath, 'analyses', 'results') + +# Create results directory if it doesn't exist +if not os.path.exists(microalgae_results_filepath): + os.makedirs(microalgae_results_filepath) + +system = microalgae_sys = microalgae_mcca_sys +tea = microalgae_tea + +# Create model with metrics +def create_model(): + """Create evaluation model with metrics for microalgae system""" + # Get system components + u = microalgae_mcca_sys.flowsheet.unit + s = microalgae_mcca_sys.flowsheet.stream + + # Find main product and boiler + main_product = s.caproic_acid_product + main_product_chemical_IDs = ['CaproicAcid'] + + # Find boiler + boiler = None + for unit in microalgae_mcca_sys.units: + if hasattr(unit, 'natural_gas') or 'BT' in unit.ID: + boiler = unit + break + + # Create LCA object + lca = create_microalgae_lca(microalgae_mcca_sys, main_product, main_product_chemical_IDs, boiler) + + # Define metrics + metrics = [ + Metric('MPSP', lambda: microalgae_tea.solve_price(main_product), '$/kg'), + Metric('TCI', lambda: microalgae_tea.TCI/1e6, 'MM$'), + Metric('VOC', lambda: microalgae_tea.VOC/1e6, 'MM$/y'), + Metric('FOC', lambda: microalgae_tea.FOC/1e6, 'MM$/y'), + Metric('GWP', lambda: lca.GWP, 'kg CO2-eq/kg'), + Metric('FEC', lambda: lca.FEC, 'MJ/kg'), + ] + + # Create namespace for parameter loading + namespace_dict = { + 'microalgae_sys': microalgae_mcca_sys, + 'microalgae_tea': microalgae_tea, + 'u': u, + 's': s, + 'lca': lca, + 'bst': bst, + 'np': np, + # Add chemical streams for easier access + 'microalgae': None, # Will be set dynamically + 'GlucoAmylase': None, # Will be set dynamically + 'AlphaAmylase': None, # Will be set dynamically + 'Yeast': None, + 'OleylAlcohol': None, + 'base_fermentation': None, + 'FGD_lime': None, + 'PowerUtility': bst.PowerUtility, + } + + # Try to find chemical streams + for stream in microalgae_mcca_sys.feeds: + if 'microalgae' in stream.ID.lower(): + namespace_dict['microalgae'] = stream + break + + # Create model + model = MicroalgaeModel(microalgae_mcca_sys, metrics=metrics, namespace_dict=namespace_dict) + + return model, lca, namespace_dict + +model, lca, namespace_dict = create_model() + +def get_adjusted_MSP(): + """Get adjusted minimum selling price""" + return microalgae_tea.solve_price(microalgae_sys.flowsheet.stream.caproic_acid_product) + +# %% + +N_simulations_per_mode = 2000 # 2000 + +percentiles = [0, 0.05, 0.1, 0.25, 0.5, 0.75, 0.9, 0.95, 1] + +notification_interval = 10 + +results_dict = {'Baseline':{'MPSP':{}, 'GWP100a':{}, 'FEC':{}, + 'GWP Breakdown':{}, 'FEC Breakdown':{},}, + 'Uncertainty':{'MPSP':{}, 'GWP100a':{}, 'FEC':{}}, + 'Sensitivity':{'Spearman':{'MPSP':{}, 'GWP100a':{}, 'FEC':{}}},} + +modes = [ + 'baseline', + ] + +parameter_distributions_filenames = [ + 'parameter_distributions.xlsx', + ] + +#%% + +timer = TicToc('timer') +timer.tic() + +# Set seed to make sure each time the same set of random numbers will be used +np.random.seed(3221) # 3221 + +for i in range(len(modes)): + mode = modes[i] + parameter_distributions_filename = os.path.join(microalgae_filepath, parameter_distributions_filenames[i]) + + print(f'\n\nLoading parameter distributions ({mode}) ...') + model.parameters = () + model.load_parameter_distributions(parameter_distributions_filename, namespace_dict) + print(f'\nLoaded parameter distributions ({mode}).') + + parameters = model.get_parameters() + + print('\n\nLoading samples ...') + samples = model.sample(N=N_simulations_per_mode, rule='L') + model.load_samples(samples) + print('\nLoaded samples.') + + model.exception_hook = 'warn' + print('\n\nSimulating baseline ...') + baseline_initial = model.metrics_at_baseline() + baseline = pd.DataFrame(data=np.array([[i for i in baseline_initial.values],]), + columns=baseline_initial.keys()) + + results_dict['Baseline']['MPSP'][mode] = get_adjusted_MSP() + results_dict['Baseline']['GWP100a'][mode] = tot_GWP = lca.GWP + results_dict['Baseline']['FEC'][mode] = tot_FEC = lca.FEC + + # GWP breakdown analysis + try: + material_GWP_breakdown = lca.material_GWP_breakdown + + results_dict['Baseline']['GWP Breakdown'][mode] = { + 'feedstock': lca.feedstock_GWP, + 'material inputs': lca.material_GWP, + 'natural gas\n(for steam generation)': getattr(lca, 'ng_GWP', 0), + 'net electricity': lca.net_electricity_GWP, + 'direct non-biogenic\nemissions': getattr(lca, 'direct_emissions_GWP', 0), + } + + tot_positive_GWP = sum([v for v in results_dict['Baseline']['GWP Breakdown'][mode].values() if v>0]) + if tot_positive_GWP > 0: + for k, v in results_dict['Baseline']['GWP Breakdown'][mode].items(): + results_dict['Baseline']['GWP Breakdown'][mode][k] = v/tot_positive_GWP + except Exception as e: + print(f"Warning: Could not calculate GWP breakdown: {e}") + results_dict['Baseline']['GWP Breakdown'][mode] = {} + + # FEC breakdown analysis + try: + material_FEC_breakdown = lca.material_FEC_breakdown + + results_dict['Baseline']['FEC Breakdown'][mode] = { + 'feedstock': lca.feedstock_FEC, + 'material inputs': lca.material_FEC, + 'natural gas\n(for steam generation)': getattr(lca, 'ng_FEC', 0), + 'net electricity': lca.net_electricity_FEC, + } + + tot_positive_FEC = sum([v for v in results_dict['Baseline']['FEC Breakdown'][mode].values() if v>0]) + if tot_positive_FEC > 0: + for k, v in results_dict['Baseline']['FEC Breakdown'][mode].items(): + results_dict['Baseline']['FEC Breakdown'][mode][k] = v/tot_positive_FEC + except Exception as e: + print(f"Warning: Could not calculate FEC breakdown: {e}") + results_dict['Baseline']['FEC Breakdown'][mode] = {} + + print(f"\nSimulated baseline. MPSP = ${round(results_dict['Baseline']['MPSP'][mode],2)}/kg.") + print('\n\nEvaluating ...') + model.evaluate(notify=notification_interval, autoload=None, autosave=None, file=None) + print('\nFinished evaluation.') + + # Baseline results + print('\n\nRe-simulating baseline ...') + baseline_end = model.metrics_at_baseline() + print(f"\nRe-simulated baseline. MPSP = ${round(results_dict['Baseline']['MPSP'][mode],2)}/kg.") + dateTimeObj = datetime.now() + minute = '0' + str(dateTimeObj.minute) if len(str(dateTimeObj.minute))==1 else str(dateTimeObj.minute) + file_to_save = os.path.join(microalgae_results_filepath, + f'_microalgae_{dateTimeObj.year}.{dateTimeObj.month}.{dateTimeObj.day}-{dateTimeObj.hour}.{minute}'\ + + f'_{N_simulations_per_mode}sims') + + baseline.index = ('initial', ) + baseline.to_excel(file_to_save+'_'+mode+'_0_baseline.xlsx') + + # Parameters + parameters = model.get_parameters() + index_parameters = len(model.get_baseline_sample()) + parameter_values = model.table.iloc[:, :index_parameters].copy() + + #%% + + # TEA results + for index_TEA, i in enumerate(model.metrics): + if hasattr(i, 'element') and i.element == 'LCA': + break + else: + index_TEA = len(model.metrics) - 2 # Assume last 2 are LCA metrics + + index_TEA = index_parameters + index_TEA + TEA_results = model.table.iloc[:, index_parameters:index_TEA].copy() + TEA_percentiles = TEA_results.quantile(q=percentiles) + + # LCA_results + LCA_results = model.table.iloc[:, index_TEA::].copy() + LCA_percentiles = LCA_results.quantile(q=percentiles) + + # # Spearman's rank correlation + table = model.table + model.table = model.table.dropna() + + spearman_results = model.spearman() + spearman_results.columns = pd.Index([i.name_with_units for i in model.metrics]) + + model.table = table + + # Calculate the cumulative probabilities of each parameter + probabilities = {} + for i in range(index_parameters): + p = parameters[i] + p_values = parameter_values.iloc[:, 2*i] + probabilities[p.name] = p.distribution.cdf(p_values) + parameter_values.insert(loc=2*i+1, + column=(parameter_values.iloc[:, 2*i].name[0], 'Probability'), + value=probabilities[p.name], + allow_duplicates=True) + + run_number = samples.shape[0] + + #%% + '''Output to Excel''' + with pd.ExcelWriter(file_to_save+'_'+mode+'_1_full_evaluation.xlsx') as writer: + parameter_values.to_excel(writer, sheet_name='Parameters') + TEA_results.to_excel(writer, sheet_name='TEA results') + TEA_percentiles.to_excel(writer, sheet_name='TEA percentiles') + LCA_results.to_excel(writer, sheet_name='LCA results') + LCA_percentiles.to_excel(writer, sheet_name='LCA percentiles') + spearman_results.to_excel(writer, sheet_name='Spearman') + model.table.to_excel(writer, sheet_name='Raw data') + + # Extract results for plotting + # Find the MPSP column + mpsp_col = None + for col in model.table.columns: + if 'MPSP' in str(col) or 'price' in str(col).lower(): + mpsp_col = col + break + + # Find GWP column + gwp_col = None + for col in model.table.columns: + if 'GWP' in str(col): + gwp_col = col + break + + # Find FEC column + fec_col = None + for col in model.table.columns: + if 'FEC' in str(col): + fec_col = col + break + + if mpsp_col is not None: + results_dict['Uncertainty']['MPSP'][mode] = model.table[mpsp_col] + else: + print("Warning: Could not find MPSP column in results") + results_dict['Uncertainty']['MPSP'][mode] = [results_dict['Baseline']['MPSP'][mode]] * len(model.table) + + if gwp_col is not None: + results_dict['Uncertainty']['GWP100a'][mode] = model.table[gwp_col] + else: + print("Warning: Could not find GWP column in results") + results_dict['Uncertainty']['GWP100a'][mode] = [results_dict['Baseline']['GWP100a'][mode]] * len(model.table) + + if fec_col is not None: + results_dict['Uncertainty']['FEC'][mode] = model.table[fec_col] + else: + print("Warning: Could not find FEC column in results") + results_dict['Uncertainty']['FEC'][mode] = [results_dict['Baseline']['FEC'][mode]] * len(model.table) + + # Spearman correlations for sensitivity analysis + df_rho, df_p = model.spearman_r() + + if mpsp_col is not None and mpsp_col in df_rho.columns: + results_dict['Sensitivity']['Spearman']['MPSP'][mode] = df_rho[mpsp_col] + else: + results_dict['Sensitivity']['Spearman']['MPSP'][mode] = pd.Series() + + if gwp_col is not None and gwp_col in df_rho.columns: + results_dict['Sensitivity']['Spearman']['GWP100a'][mode] = df_rho[gwp_col] + else: + results_dict['Sensitivity']['Spearman']['GWP100a'][mode] = pd.Series() + + if fec_col is not None and fec_col in df_rho.columns: + results_dict['Sensitivity']['Spearman']['FEC'][mode] = df_rho[fec_col] + else: + results_dict['Sensitivity']['Spearman']['FEC'][mode] = pd.Series() + +#%% Clean up NaN values for plotting +metrics = ['MPSP', 'GWP100a', 'FEC'] +tot_NaN_vals_dict = results_dict['Errors'] = {metric: {mode: 0 for mode in modes} for metric in metrics} +for mode in modes: + for metric in metrics: + median_val = 1.5 # Default fallback value + if len(results_dict['Uncertainty'][metric][mode]) > 0: + median_val = np.nanmedian(results_dict['Uncertainty'][metric][mode]) + if np.isnan(median_val): + median_val = 1.5 + + for i in range(len(results_dict['Uncertainty'][metric][mode])): + if np.isnan(results_dict['Uncertainty'][metric][mode].iloc[i]): + results_dict['Uncertainty'][metric][mode].iloc[i] = median_val + tot_NaN_vals_dict[metric][mode] += 1 + +# %% Plots - temporarily disabled due to contourplots dependency +# MPSP_units = r"$\mathrm{\$}\cdot\mathrm{kg}^{-1}$" +# GWP_units = r"$\mathrm{kg}$"+" "+ r"$\mathrm{CO}_{2}\mathrm{-eq.}\cdot\mathrm{kg}^{-1}$" +# FEC_units = r"$\mathrm{MJ}\cdot\mathrm{kg}^{-1}$" + +scenario_name_labels = ['Baseline'] + +def get_small_range(num, offset): + return(num-offset, num+offset) + +print(f'\nAnalysis completed. Timer: {timer.toc():.2f} seconds') +print(f'Results saved to: {file_to_save}') + +# Print summary +print('\n=== UNCERTAINTY ANALYSIS SUMMARY ===') +for mode in modes: + print(f'\n{mode.upper()} MODE:') + print(f' MPSP: {results_dict["Baseline"]["MPSP"][mode]:.3f} $/kg') + print(f' GWP: {results_dict["Baseline"]["GWP100a"][mode]:.3f} kg CO2-eq/kg') + print(f' FEC: {results_dict["Baseline"]["FEC"][mode]:.3f} MJ/kg') + + if len(results_dict['Uncertainty']['MPSP'][mode]) > 0: + print(f' MPSP range: {np.min(results_dict["Uncertainty"]["MPSP"][mode]):.3f} - {np.max(results_dict["Uncertainty"]["MPSP"][mode]):.3f} $/kg') + print(f' GWP range: {np.min(results_dict["Uncertainty"]["GWP100a"][mode]):.3f} - {np.max(results_dict["Uncertainty"]["GWP100a"][mode]):.3f} kg CO2-eq/kg') + print(f' FEC range: {np.min(results_dict["Uncertainty"]["FEC"][mode]):.3f} - {np.max(results_dict["Uncertainty"]["FEC"][mode]):.3f} MJ/kg') + +# Print detailed breakdown analysis +print('\n=== DETAILED BREAKDOWN ANALYSIS ===') +for mode in modes: + print(f'\n{mode.upper()} MODE BREAKDOWN:') + + # GWP Breakdown + if results_dict['Baseline']['GWP Breakdown'][mode]: + print(f'\nGWP Breakdown (Total: {results_dict["Baseline"]["GWP100a"][mode]:.3f} kg CO2-eq/kg):') + gwp_breakdown = results_dict['Baseline']['GWP Breakdown'][mode] + total_gwp = sum([abs(v) for v in gwp_breakdown.values()]) + for component, value in gwp_breakdown.items(): + percentage = (value / total_gwp * 100) if total_gwp > 0 else 0 + print(f' {component}: {value:.4f} kg CO2-eq/kg ({percentage:.1f}%)') + + # FEC Breakdown + if results_dict['Baseline']['FEC Breakdown'][mode]: + print(f'\nFEC Breakdown (Total: {results_dict["Baseline"]["FEC"][mode]:.3f} MJ/kg):') + fec_breakdown = results_dict['Baseline']['FEC Breakdown'][mode] + total_fec = sum([abs(v) for v in fec_breakdown.values()]) + for component, value in fec_breakdown.items(): + percentage = (value / total_fec * 100) if total_fec > 0 else 0 + print(f' {component}: {value:.4f} MJ/kg ({percentage:.1f}%)') + +# Print sensitivity analysis results +print('\n=== SENSITIVITY ANALYSIS (SPEARMAN CORRELATIONS) ===') +for mode in modes: + print(f'\n{mode.upper()} MODE SENSITIVITY:') + + # MPSP Correlations + if not results_dict['Sensitivity']['Spearman']['MPSP'][mode].empty: + print(f'\nTop 10 correlations with MPSP:') + mpsp_corr = results_dict['Sensitivity']['Spearman']['MPSP'][mode].copy() + # Remove NaN values and sort by absolute correlation + mpsp_corr = mpsp_corr.dropna() + if len(mpsp_corr) > 0: + sorted_corr = mpsp_corr.abs().sort_values(ascending=False) + for i, (param, abs_corr) in enumerate(sorted_corr.head(10).items()): + actual_corr = mpsp_corr[param] + print(f' {i+1:2d}. {param}: {actual_corr:.3f}') + else: + print(' No significant correlations found') + + # GWP Correlations + if not results_dict['Sensitivity']['Spearman']['GWP100a'][mode].empty: + print(f'\nTop 10 correlations with GWP:') + gwp_corr = results_dict['Sensitivity']['Spearman']['GWP100a'][mode].copy() + gwp_corr = gwp_corr.dropna() + if len(gwp_corr) > 0: + sorted_corr = gwp_corr.abs().sort_values(ascending=False) + for i, (param, abs_corr) in enumerate(sorted_corr.head(10).items()): + actual_corr = gwp_corr[param] + print(f' {i+1:2d}. {param}: {actual_corr:.3f}') + else: + print(' No significant correlations found') + + # FEC Correlations + if not results_dict['Sensitivity']['Spearman']['FEC'][mode].empty: + print(f'\nTop 10 correlations with FEC:') + fec_corr = results_dict['Sensitivity']['Spearman']['FEC'][mode].copy() + fec_corr = fec_corr.dropna() + if len(fec_corr) > 0: + sorted_corr = fec_corr.abs().sort_values(ascending=False) + for i, (param, abs_corr) in enumerate(sorted_corr.head(10).items()): + actual_corr = fec_corr[param] + print(f' {i+1:2d}. {param}: {actual_corr:.3f}') + else: + print(' No significant correlations found') + +# Print statistics summary +print('\n=== STATISTICAL SUMMARY ===') +for mode in modes: + print(f'\n{mode.upper()} MODE STATISTICS:') + + for metric in ['MPSP', 'GWP100a', 'FEC']: + if len(results_dict['Uncertainty'][metric][mode]) > 0: + data = results_dict['Uncertainty'][metric][mode] + mean_val = np.mean(data) + std_val = np.std(data) + p5 = np.percentile(data, 5) + p95 = np.percentile(data, 95) + median_val = np.median(data) + + units = {'MPSP': '$/kg', 'GWP100a': 'kg CO2-eq/kg', 'FEC': 'MJ/kg'} + unit = units[metric] + + print(f'\n{metric}:') + print(f' Mean: {mean_val:.3f} {unit}') + print(f' Std Dev: {std_val:.3f} {unit}') + print(f' Median: {median_val:.3f} {unit}') + print(f' 5th percentile: {p5:.3f} {unit}') + print(f' 95th percentile: {p95:.3f} {unit}') + print(f' Range: {np.min(data):.3f} - {np.max(data):.3f} {unit}') + +# Print error summary +print('\n=== ERROR SUMMARY ===') +total_errors = 0 +for metric in metrics: + for mode in modes: + errors = tot_NaN_vals_dict[metric][mode] + if errors > 0: + print(f'{metric} ({mode}): {errors} NaN values replaced') + total_errors += errors + +if total_errors == 0: + print('No errors encountered during simulation.') +else: + print(f'Total errors handled: {total_errors}') + +print(f'\n=== ANALYSIS COMPLETED ===') +print(f'Total simulation time: {timer.toc():.2f} seconds') +print(f'Results saved to: {file_to_save}') +print(f'Number of parameters analyzed: {len(parameters)}') +print(f'Number of successful simulations: {len(model.table)}') +print(f'Output files generated:') +print(f' - Baseline: {file_to_save}_{mode}_0_baseline.xlsx') +print(f' - Full results: {file_to_save}_{mode}_1_full_evaluation.xlsx') + +if __name__ == '__main__': + print("Microalgae uncertainty analysis completed!") \ No newline at end of file diff --git a/biorefineries/microalgae/units.py b/biorefineries/microalgae/units.py new file mode 100644 index 00000000..6b2d4003 --- /dev/null +++ b/biorefineries/microalgae/units.py @@ -0,0 +1,831 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Sat June 23 17:41:00 2025 + +Microalgae biorefinery to produce medium chain fatty acids +by anaerobic fermentation without external electron donor addition- chemicals database + +References +---------- +[1] BioSTEAM Documentation: + https://biosteam.readthedocs.io/en/latest/API/thermosteam/Chemicals.html +[2] Cortes-Peña et al., BioSTEAM: A Fast and Flexible Platform for the Design, + Simulation, and Techno-Economic Analysis of Biorefineries under Uncertainty. + ACS Sustainable Chem. Eng. 2020, 8 (8), 3302–3310. +[3] 3-Hydroxypropionic acid biorefineries project: + https://github.com/BioSTEAMDevelopmentGroup/Bioindustrial-Park/tree/master/biorefineries/HP + +@author: Xingdong Shi +@version: 0.0.1 +""" + +import numpy as np +import biosteam as bst +import math +from biosteam import Unit +from biosteam.units.design_tools import pressure_vessel_material_factors +from biosteam.units.decorators import cost +from biosteam.units import HXutility, Mixer, MixTank, Pump, StorageTank, StirredTankReactor +from thermosteam import Stream +import thermosteam as tmo +from ._chemicals import chems +from .utils import baseline_feedflow, CEPCI +from thermosteam import reaction as rxn +from biorefineries.lactic._units import Reactor + +def anaerobic_rxn(reactant): + MW = getattr(chems, reactant).MW + return rxn.Reaction(f"{1/MW}{reactant} -> {P_ch4}CH4 + {P_co2}CO2 + {P_sludge}WWTsludge", reactant, 0.91) + +def find_split(*tuples): + + split = {} + for ID, flow0, flow1 in tuples: + total = flow0 + flow1 + split[ID] = flow0 / total if total else 0 + return split + +Rxn = tmo.reaction.Reaction +ParallelRxn = tmo.reaction.ParallelReaction +ln = math.log +exp = math.exp +_gal2m3 = 0.003785 +_gpm2m3hr = 0.227124 +chems = tmo.settings.chemicals +organics = ['Glucose', 'Xylose', 'Protein', 'Ash', 'AceticAcid', 'Ethanol'] +if 'WWTsludge' in organics: + organics.remove('WWTsludge') +P_sludge = 0.05/0.91/chems.WWTsludge.MW +MW = np.array([chems.CH4.MW, chems.CO2.MW]) +CH4_molcomp = 0.60 +mass = np.array([CH4_molcomp, 1-CH4_molcomp])*MW +mass /= mass.sum() +mass *= 0.381/(0.91) +P_ch4, P_co2 = mass/MW +anaerobic_digestion = rxn.ParallelReaction([anaerobic_rxn(i) for i in organics]) + +splits = [ + ('Glucose', 3, 42), + ('Xylose', 7, 85), + ('Protein', 51, 18), + ('Ash', 813, 280), + ('AceticAcid', 0, 5), + ('Ethanol', 1, 15), +] + +tmo.settings.set_thermo(chems) + +# %% +# ========================= +# Microalgae Crushing +# ========================= +@cost(basis='Flow rate', ID='System', units='kg/hr', + kW=511.3205, cost=13329690, S=94697, CE=CEPCI[2009], n=0.6, BM=1.7) +class FeedstockPreprocessing(Unit): + _baseline_flow_rate = baseline_feedflow.sum() + _cached_flow_rate = 2205 + _N_ins = 1 + _N_outs = 1 + + def _run(self): + self.outs[0].copy_like(self.ins[0]) + s = self.outs[0] + for comp in s.chemicals.IDs: + if comp != 'Microalgae': + s.imass[comp] = 0 + +# %% +# ============================================================================= +# Pretreatment +# ============================================================================= +@cost(basis='Flow rate', ID='Tank', units='kg/hr', + cost=6210, S=1981, CE=CEPCI[2010], n=0.7, BM=3) +@cost(basis='Flow rate', ID='Pump', units='kg/hr', + cost=8000, S=3720, CE=CEPCI[2009], n=0.8, BM=2.3) +class SulfuricAcidAdditionTank(StorageTank): + _N_ins = 1 + _N_outs = 1 + # acid_volume_per_biomass = 0.8 mL/g + # acid_density = 1.84 # g/mL + acid_loading = 1.47 # g H2SO4 / g microalgae + + def _run(self): + feed = self.ins[0] + out = self.outs[0] + microalgae_mass = feed.imass['Microalgae'] + + pure_H2SO4 = microalgae_mass * self.acid_loading + + acid_purity = 0.93 + acid_solution_mass = pure_H2SO4 / acid_purity + water_mass = acid_solution_mass * (1 - acid_purity) + out.copy_like(feed) + out.imass['H2SO4'] = pure_H2SO4 + out.imass['H2O'] = water_mass + +# Sulfuric acid in-line mixer +@cost(basis='Flow rate', ID='Mixer', units='kg/hr', + cost=6000, S=136260, CE=CEPCI[2009], n=0.5, BM=1) +class SulfuricAcidMixer(Unit): + _N_ins = 2 + _N_outs = 1 + def _run(self): + microalgae, acid = self.ins + out = self.outs[0] + out.mix_from([microalgae, acid]) + +# Pretreatment reactor +@cost(basis='Dry flow rate', ID='Pretreatment reactor', units='kg/hr', + kW=5120, cost=19812400, S=83333, CE=CEPCI[2009], n=0.6, BM=1.5) +class AcidPretreatmentReactor(Unit): + _N_ins = 1 + _N_outs = 1 + T = 121 + 273.15 # K + t = 20 / 60 # h + + microalgae_decomposition = tmo.reaction.Reaction( + 'Microalgae -> 0.5 Protein + 0.159 Lipid + 0.234 Carbohydrate + 0.107 Ash', + reactant='Microalgae', + X=1.0, + basis='wt' + ) + + def _run(self): + feed = self.ins[0] + out = self.outs[0] + out.copy_like(feed) + self.microalgae_decomposition(out) + +# pH adjustment +@cost(basis='Flow rate', ID='Mixer', units='kg/hr', + # Size basis on the total flow, not just ammonia, + # thus assuming difference caused by MWs of NH3 and NH4OH is negligible + cost=50000, S=157478, CE=CEPCI[2009], n=0.5, BM=1.5) +class NAOHMixer(Mixer): + def _run(self): + NAOH, water = self.ins + mixture = self.outs[0] + + # Make 10% solution + water.imass['Water'] = NAOH.imass['NAOH'] * (1-0.1)/0.1 + mixture.mix_from([NAOH, water]) + + +@cost(basis='Flow rate', ID='Tank', units='kg/hr', + cost=236000, S=410369, CE=CEPCI[2009], n=0.7, BM=2) +@cost(basis='Flow rate', ID='Agitator', units='kg/hr', + kW=7.457, cost=21900, S=410369, CE=CEPCI[2009], n=0.5, BM=1.5) +class NAOHAdditionTank(MixTank): + # Based on experimental data: 0.64kg NaOH consumed per ton of microalgae + naoh_mass_per_ton_microalgae = 0.64 # kg NaOH per ton microalgae + + def _run(self): + # Calculate required NaOH based on microalgae mass + if hasattr(self, 'biomass_mass') and self.biomass_mass > 0: + required_naoh = self.biomass_mass * self.naoh_mass_per_ton_microalgae / 1000 # convert to kg/hr + # Set NaOH input amount + if len(self.ins) > 0: + self.ins[0].imass['NaOH'] = required_naoh + + # Reaction definition Reactant Conversion + self.neutralization_rxn = Rxn('2 NaOH + H2SO4 -> Na2SO4 + 2 H2O', 'H2SO4', 1) + + MixTank._run(self) + self.neutralization_rxn.adiabatic_reaction(self.outs[0]) + +@cost(basis='Flow rate', ID='Pump', units='kg/hr', + kW=74.57, cost=22500, S=402194, CE=CEPCI[2009], n=0.8, BM=2.3) +class HydrolysatePump(Unit): + _units = {'Flow rate': 'kg/hr'} + _graphics = Pump._graphics + + def _design(self): + Design = self.design_results + Design['Flow rate'] = self.outs[0].F_mass + + + +# %% +# ============================================================================= +# Conversion +# ============================================================================= +# Enzyme hydrolysate mixer +@cost(basis='Flow rate', ID='Mixer', units='kg/hr', + kW=74.57, cost=109000, S=379938, CE=CEPCI[2009], n=0.5, BM=1.7) +class EnzymeHydrolysateMixer(Mixer): + _N_ins = 3 + _N_outs = 1 + + # Experimental data: 6kg GlucoAmylase and 266.7kg AlphaAmylase consumed per ton of microalgae + glucoamylase_per_ton = 6 # kg per ton microalgae + alphaamylase_per_ton = 266.7 # kg per ton microalgae + # enzyme_mass_per_ton_microalgae = 172.7 # kg enzyme per ton microalgae (deprecated) + + def _run(self): + hydrolysate, enzyme, water = self.ins + effluent = self.outs[0] + + # Calculate consumption of two enzymes based on microalgae mass + if hasattr(self, 'biomass_mass') and self.biomass_mass > 0: + required_glucoamylase = self.biomass_mass * self.glucoamylase_per_ton / 1000 # kg/hr + required_alphaamylase = self.biomass_mass * self.alphaamylase_per_ton / 1000 # kg/hr + enzyme.imass['GlucoAmylase'] = required_glucoamylase + enzyme.imass['AlphaAmylase'] = required_alphaamylase + else: + # If no biomass_mass attribute, use default method + enzyme.imass['GlucoAmylase'] = hydrolysate.F_mass * 0.006 # 0.6% + enzyme.imass['AlphaAmylase'] = hydrolysate.F_mass * 0.2667 # 26.67% + + effluent.mix_from([hydrolysate, enzyme, water]) + + + +@cost(basis='Saccharification tank size', ID='Saccharification tanks', units='kg', + cost=3840000, S=421776*24, CE=CEPCI[2009], n=0.7, BM=2) +@cost(basis='Slurry flow rate', ID='Saccharification transfer pumps', units='kg/hr', + kW=74.57, cost=47200, S=421776*24, CE=CEPCI[2009], n=0.8, BM=2.3) +class Saccharification(Unit): + _N_ins = 1 + _N_outs = 1 + _units= {'Saccharification tank size': 'kg', + 'Slurry flow rate': 'kg/hr'} + + # Residence time of countinuous saccharification tanks (hr) + tau_saccharification = 3 + + def _init(self, ID='', ins=None, outs=(), T=50+273.15): + + self.T = T + + # Based on Table 9 on Page 28 of Humbird et al. + #!!! Yalin updated the glucan-to-glucose conversion to match the table + self.saccharification_rxns = ParallelRxn([ + # Reaction definition Reactant Conversion + Rxn('Carbohydrate + H2O -> Glucose', 'Carbohydrate', 1) + ]) + + # self.saccharified_stream = Stream(None) + + def _run(self): + feed = self.ins[0] + ss= self.outs[0] + ss.copy_like(feed) + ss.T = self.T + self.saccharification_rxns(ss.mol) + neutralization_rxn = Rxn('2 NaOH + H2SO4 -> Na2SO4 + 2 H2O', 'NaOH', 1) + neutralization_rxn(ss) + + + + def _design(self): + Design = self.design_results + total_mass_flow = self.ins[0].F_mass + Design['Saccharification tank size'] = total_mass_flow * self.tau_saccharification + Design['Slurry flow rate'] = total_mass_flow + +## MCCA Fermentation +@cost(basis='Fermenter size', ID='Fermenter', units='kg', + kW=100, cost=5128000, S=(42607+443391+948+116)*(60+36), + CE=CEPCI[2009], n=0.7, BM=1.5) +@cost(basis='Fermenter size', ID='Fermenter agitator', units='kg', + kW=10, + cost=630000, + S=(42607+443391+948+116)*(60+36), + CE=CEPCI[2009], n=0.7, BM=1.5) + +class MCCAFermentation(StirredTankReactor): + _N_ins = 5 + _N_outs = 4 + + auxiliary_unit_names = ('heat_exchanger',) + _units= { + **Reactor._units, + 'Fermenter size': 'kg', + 'Broth flow rate': 'kg/hr', + } + + tau_batch_turnaround = 12 # in hr, seedmentaiton time + tau_cofermentation = 15 * 24 # in hr, 15 d + mol_NH4OH_per_acid_pH_control = 0.2988 # from IBRL batch run + mol_lime_per_acid_pH_control = 0.2988/2 # mol NH4OH / 2 + pH = 5.0 + + def __init__(self, ID='', ins=None, outs=(), thermo=None, *, T=37+273.15, + P=101325, V_wf=0.8, length_to_diameter=0.6, + kW_per_m3=0.02 * 0.7457 / 0.075, # ref from succinic projects + mixing_intensity=300, + #wall_thickness_factor=1, + vessel_material='Stainless steel 304', + vessel_type='Vertical', + neutralization=True, + neutralizing_agent='Lime', + mode='batch', feed_freq=1, + pH_control=True, + base_neutralizes_product=True, + tau = tau_cofermentation, + microalgae_mass_flow: float | None = None, # kg/hr basis of original algae for yield calc + titer: float = 2.003, # g/L + # allow_dilution=False, + # allow_concentration=False, + # sugars=None + ): + + StirredTankReactor.__init__(self, ID, ins, outs) + + self.T = T + self.P = P + self.V_wf = V_wf + self.length_to_diameter = length_to_diameter + self.mixing_intensity = mixing_intensity + self.kW_per_m3 = kW_per_m3 + #self.wall_thickness_factor = wall_thickness_factor + self.vessel_material = vessel_material + self.vessel_type = vessel_type + self.neutralization = neutralization + self.mode = mode + self.feed_freq = feed_freq + self.pH_control=pH_control + self.base_neutralizes_product = base_neutralizes_product + self.neutralization = neutralization + self.neutralizing_agent = neutralizing_agent + self.tau = tau + + #self.allow_dilution = allow_dilution + #self.allow_concentration = allow_concentration + #self.sugars = sugars or tuple(i.ID for i in self.chemicals.sugars) + + # Store optional explicit biomass flow for yield calculations + self.microalgae_mass_flow = microalgae_mass_flow + self.titer = titer + + ID = self.ID + self._mixed_feed = Stream(f'{ID}_mixed_feed') + self._tot_feed = Stream(f'{ID}_tot_feed') + # Before reaction, after reaction, with last feed + self._single_feed0 = Stream(f'{ID}_single_feed0') + self._single_feed1 = Stream(f'{ID}_single_feed1') + self._last = Stream(f'{ID}_last') + self._init = Stream(f'{ID}_init') + hx_in = Stream(f'{ID}_hx_in') + hx_out = Stream(f'{ID}_hx_out') + # Add '.' in ID for auxiliary units + self.heat_exchanger = HXutility(ID=f'.{ID}_hx', ins=hx_in, outs=hx_out, T=T) + + # the reaction is simplified for simulation + #self.mcca_rxns = ParallelRxn([ + #Rxn('Protein -> AceticAcid', 'Protein', 0.1), + #Rxn('Protein -> PropionicAcid', 'Protein', 0.1), + #Rxn('Protein -> ButyricAcid', 'Protein', 0.1), + #Rxn('Protein -> CaproicAcid', 'Protein', 0.1), + #Rxn('Glucose -> Ethanol', 'Glucose', 0.018), + #Rxn('Glucose -> Butanol', 'Glucose', 0.0022), + #Rxn('Glucose -> AceticAcid', 'Glucose', 0.10094), + #Rxn('Glucose -> PropionicAcid', 'Glucose', 0.00977), + #Rxn('Glucose -> ValericAcid', 'Glucose', 0.012), + #Rxn('Glucose -> HeptanoicAcid', 'Glucose', 0.004), + #Rxn('Glucose -> ButyricAcid', 'Glucose', 0.122), + #Rxn('Glucose -> CaproicAcid', 'Glucose', 0.134), + #Rxn('Glucose -> CaprylicAcid', 'Glucose', 0.008), + #]) + + # This is a Siplified Neutralization, because we don't want to affect the products yield + self.lime_neutralization_rxns = ParallelRxn([ + # Reaction definition Reactant Conversion + Rxn('H2SO4 + Lime -> CaSO4 + 2 H2O', 'H2SO4', 1) + ]) + + self.lime_pH_control_rxns = ParallelRxn([ + # Reaction definition Reactant Conversion + Rxn('H2SO4 + Lime -> CaSO4 + 2 H2O', 'H2SO4', 1) + ]) + + + # self.mcca_rxns = ParallelRxn([ + # Rxn('Glucose -> 2 Ethanol + 2 CO2', 'Glucose', 0.7), + # Rxn('Glucose -> 2 AceticAcid + 2 CO2 + 2 H2', 'Glucose', 0.2), + # Rxn('Glucose -> 2 PropionicAcid + 2 CO2 + 2 H2', 'Glucose', 0.1), + # Rxn('5 PropionicAcid + 6 Ethanol -> 5 ValericAcid + AceticAcid + 4 H2O + 2 H2', 'PropionicAcid', 0.7), + # Rxn('5 ValericAcid + 6 Ethanol -> 5 HeptanoicAcid + AceticAcid + 4 H2O + 2 H2', 'ValericAcid', 0.7), + # Rxn('ButyricAcid + Ethanol -> AceticAcid + Butanol', 'ButyricAcid', 0.1), + # Rxn('4 AceticAcid + 6 Ethanol -> 5 ButyricAcid + 4 H2O + 2 H2', 'AceticAcid', 0.7), + # Rxn('5 ButyricAcid + 6 Ethanol -> 5 CaproicAcid + AceticAcid + 4 H2O + 2 H2', 'ButyricAcid', 0.7), + # Rxn('5 CaproicAcid + 6 Ethanol -> 5 CaprylicAcid + AceticAcid + 4 H2O + 2 H2', 'CaproicAcid', 0.2), + # ]) + + def _run(self): + substrate, yeast_seed, FermMicrobe, lime, n2 = self.ins + broth, gas, fermentation_waste, yeast_recycle = self.outs + broth.mix_from(self.ins) + broth.T = self.T + broth.P = self.P + # Determine basis biomass mass (kg/hr) + microalgae_mass = self.microalgae_mass_flow + caproic_acid_yield_factor = getattr(self, 'caproic_acid_yield_factor', 1.0) + titer = getattr(self, 'titer', 2.003) # g/L + + base_titer = 2.003 + inhibition_titer = 10.0 + + if titer > base_titer: + if titer <= inhibition_titer: + inhibition_factor = 1.0 - 0.05 * (titer - base_titer) / (inhibition_titer - base_titer) + else: + inhibition_factor = 0.95 * np.exp(-0.1 * (titer - inhibition_titer)) + else: + inhibition_factor = 1.0 + 0.01 * (base_titer - titer) + + effective_yield_factor = caproic_acid_yield_factor * inhibition_factor + + base_yields = { + 'Ethanol': 0.01, + 'Butanol': 0.004, + 'AceticAcid': 0.1, + 'PropionicAcid': 0.009, + 'ButyricAcid': 0.18, + 'ValericAcid': 0.01, + 'CaproicAcid': 0.27, + 'HeptanoicAcid': 0.006, + 'CaprylicAcid': 0.04 + } + + total_base_yield = sum(base_yields.values()) + base_caproic_yield = base_yields['CaproicAcid'] + + new_caproic_yield = base_caproic_yield * effective_yield_factor + delta_caproic = new_caproic_yield - base_caproic_yield + total_other_base_yield = total_base_yield - base_caproic_yield + + target_caproic_mass = microalgae_mass * new_caproic_yield + target_volume_L = target_caproic_mass * 1000 / titer + + if abs(delta_caproic) > 1e-6 and total_other_base_yield > 1e-6: + adjustment_factor = max(0.1, (total_other_base_yield - delta_caproic) / total_other_base_yield) + for compound, base_yield in base_yields.items(): + if compound == 'CaproicAcid': + broth.imass[compound] = target_caproic_mass + else: + broth.imass[compound] = microalgae_mass * base_yield * adjustment_factor + else: + for compound, base_yield in base_yields.items(): + if compound == 'CaproicAcid': + broth.imass[compound] = target_caproic_mass + else: + broth.imass[compound] = microalgae_mass * base_yield + + current_volume = broth.F_vol # L + if current_volume > 0 and titer >= 5.0: + water_adjustment = max(0, target_volume_L - current_volume) + max_water = microalgae_mass * 10 + water_adjustment = min(water_adjustment, max_water) + broth.imass['Water'] += water_adjustment + gas.copy_like(n2) + gas.imass['H2'] = microalgae_mass * 0.1 + + # yeast recycle for simulating membrane bioreactor + yeast_recovery_rate = 0.95 + total_yeast_in_broth = broth.imass['Yeast'] + if total_yeast_in_broth > 0: + recoverable_yeast = total_yeast_in_broth * yeast_recovery_rate + yeast_recycle.empty() + yeast_recycle.imass['Yeast'] = recoverable_yeast + yeast_recycle.T = self.T + yeast_recycle.P = self.P + broth.imass['Yeast'] = total_yeast_in_broth - recoverable_yeast + else: + yeast_recycle.empty() + + fermentation_waste.empty() + + def _design(self): + Design = self.design_results + Design['Fermenter size'] = self.outs[0].F_mass * self.tau + Design['Broth flow rate'] = self.outs[0].F_mass + duty = 50 * self.F_mass_in # Heat duty + self.add_heat_utility(duty, self.T) + + +@cost(basis='Fermenter size', ID='Fermenter', units='kg', + kW=100, cost=5128000, S=(42607+443391+948+116)*(60+36), + CE=CEPCI[2009], n=0.7, BM=1.5) +@cost(basis='Fermenter size', ID='Fermenter agitator', units='kg', + kW=10, + cost=630000, + S=(42607+443391+948+116)*(60+36), + CE=CEPCI[2009], n=0.7, BM=1.5) +class MCCAFermentation_no_yeast(StirredTankReactor): + _N_ins = 4 + _N_outs = 3 + + auxiliary_unit_names = ('heat_exchanger',) + _units= { + **Reactor._units, + 'Fermenter size': 'kg', + 'Broth flow rate': 'kg/hr', + } + + tau_batch_turnaround = 12 # in hr, seedmentaiton time + tau_cofermentation = 15 * 24 # in hr, 15 d + mol_NH4OH_per_acid_pH_control = 0.2988 # from IBRL batch run + mol_lime_per_acid_pH_control = 0.2988/2 # mol NH4OH / 2 + pH = 5.0 + + def __init__(self, ID='', ins=None, outs=(), thermo=None, *, T=37+273.15, + P=101325, V_wf=0.8, length_to_diameter=0.6, + kW_per_m3=0.02 * 0.7457 / 0.075, # ref from succinic projects + mixing_intensity=300, + #wall_thickness_factor=1, + vessel_material='Stainless steel 304', + vessel_type='Vertical', + neutralization=True, + neutralizing_agent='Lime', + mode='batch', feed_freq=1, + pH_control=True, + base_neutralizes_product=True, + tau = tau_cofermentation, + microalgae_mass_flow: float | None = None, # kg/hr basis of original algae for yield calc + titer: float = 1.208, # g/L + # allow_dilution=False, + # allow_concentration=False, + # sugars=None + ): + + StirredTankReactor.__init__(self, ID, ins, outs) + + self.T = T + self.P = P + self.V_wf = V_wf + self.length_to_diameter = length_to_diameter + self.mixing_intensity = mixing_intensity + self.kW_per_m3 = kW_per_m3 + #self.wall_thickness_factor = wall_thickness_factor + self.vessel_material = vessel_material + self.vessel_type = vessel_type + self.neutralization = neutralization + self.mode = mode + self.feed_freq = feed_freq + self.pH_control=pH_control + self.base_neutralizes_product = base_neutralizes_product + self.neutralization = neutralization + self.neutralizing_agent = neutralizing_agent + self.tau = tau + + #self.allow_dilution = allow_dilution + #self.allow_concentration = allow_concentration + #self.sugars = sugars or tuple(i.ID for i in self.chemicals.sugars) + + # Store optional explicit biomass flow for yield calculations + self.microalgae_mass_flow = microalgae_mass_flow + self.titer = titer + + ID = self.ID + self._mixed_feed = Stream(f'{ID}_mixed_feed') + self._tot_feed = Stream(f'{ID}_tot_feed') + # Before reaction, after reaction, with last feed + self._single_feed0 = Stream(f'{ID}_single_feed0') + self._single_feed1 = Stream(f'{ID}_single_feed1') + self._last = Stream(f'{ID}_last') + self._init = Stream(f'{ID}_init') + hx_in = Stream(f'{ID}_hx_in') + hx_out = Stream(f'{ID}_hx_out') + # Add '.' in ID for auxiliary units + self.heat_exchanger = HXutility(ID=f'.{ID}_hx', ins=hx_in, outs=hx_out, T=T) + + # the reaction is simplified for simulation + #self.mcca_rxns = ParallelRxn([ + #Rxn('Protein -> AceticAcid', 'Protein', 0.1), + #Rxn('Protein -> PropionicAcid', 'Protein', 0.1), + #Rxn('Protein -> ButyricAcid', 'Protein', 0.1), + #Rxn('Protein -> CaproicAcid', 'Protein', 0.1), + #Rxn('Glucose -> Ethanol', 'Glucose', 0.018), + #Rxn('Glucose -> Butanol', 'Glucose', 0.0022), + #Rxn('Glucose -> AceticAcid', 'Glucose', 0.10094), + #Rxn('Glucose -> PropionicAcid', 'Glucose', 0.00977), + #Rxn('Glucose -> ValericAcid', 'Glucose', 0.012), + #Rxn('Glucose -> HeptanoicAcid', 'Glucose', 0.004), + #Rxn('Glucose -> ButyricAcid', 'Glucose', 0.122), + #Rxn('Glucose -> CaproicAcid', 'Glucose', 0.134), + #Rxn('Glucose -> CaprylicAcid', 'Glucose', 0.008), + #]) + + # This is a Siplified Neutralization, because we don't want to affect the products yield + self.lime_neutralization_rxns = ParallelRxn([ + # Reaction definition Reactant Conversion + Rxn('H2SO4 + Lime -> CaSO4 + 2 H2O', 'H2SO4', 1) + ]) + + self.lime_pH_control_rxns = ParallelRxn([ + # Reaction definition Reactant Conversion + Rxn('H2SO4 + Lime -> CaSO4 + 2 H2O', 'H2SO4', 1) + ]) + + + # self.mcca_rxns = ParallelRxn([ + # Rxn('Glucose -> 2 Ethanol + 2 CO2', 'Glucose', 0.7), + # Rxn('Glucose -> 2 AceticAcid + 2 CO2 + 2 H2', 'Glucose', 0.2), + # Rxn('Glucose -> 2 PropionicAcid + 2 CO2 + 2 H2', 'Glucose', 0.1), + # Rxn('5 PropionicAcid + 6 Ethanol -> 5 ValericAcid + AceticAcid + 4 H2O + 2 H2', 'PropionicAcid', 0.7), + # Rxn('5 ValericAcid + 6 Ethanol -> 5 HeptanoicAcid + AceticAcid + 4 H2O + 2 H2', 'ValericAcid', 0.7), + # Rxn('ButyricAcid + Ethanol -> AceticAcid + Butanol', 'ButyricAcid', 0.1), + # Rxn('4 AceticAcid + 6 Ethanol -> 5 ButyricAcid + 4 H2O + 2 H2', 'AceticAcid', 0.7), + # Rxn('5 ButyricAcid + 6 Ethanol -> 5 CaproicAcid + AceticAcid + 4 H2O + 2 H2', 'ButyricAcid', 0.7), + # Rxn('5 CaproicAcid + 6 Ethanol -> 5 CaprylicAcid + AceticAcid + 4 H2O + 2 H2', 'CaproicAcid', 0.2), + # ]) + + def _run(self): + substrate, FermMicrobe, lime, n2 = self.ins + broth, gas, fermentation_waste= self.outs + broth.mix_from(self.ins) + broth.T = self.T + broth.P = self.P + # Determine basis biomass mass (kg/hr) + microalgae_mass = self.microalgae_mass_flow + caproic_acid_yield_factor = getattr(self, 'caproic_acid_yield_factor', 1.0) + titer = getattr(self, 'titer', 1.208) # g/L + + base_titer = 1.208 + inhibition_titer = 10.0 + + if titer > base_titer: + if titer <= inhibition_titer: + inhibition_factor = 1.0 - 0.05 * (titer - base_titer) / (inhibition_titer - base_titer) + else: + inhibition_factor = 0.95 * np.exp(-0.1 * (titer - inhibition_titer)) + else: + inhibition_factor = 1.0 + 0.01 * (base_titer - titer) + + effective_yield_factor = caproic_acid_yield_factor * inhibition_factor + + base_yields = { + 'Ethanol': 0.01, + 'Butanol': 0.002, + 'AceticAcid': 0.1, + 'PropionicAcid': 0.009, + 'ButyricAcid': 0.086, + 'ValericAcid': 0.01, + 'CaproicAcid': 0.08, + 'HeptanoicAcid': 0.004, + 'CaprylicAcid': 0.005 + } + + total_base_yield = sum(base_yields.values()) + base_caproic_yield = base_yields['CaproicAcid'] + + new_caproic_yield = base_caproic_yield * effective_yield_factor + delta_caproic = new_caproic_yield - base_caproic_yield + total_other_base_yield = total_base_yield - base_caproic_yield + + target_caproic_mass = microalgae_mass * new_caproic_yield + target_volume_L = target_caproic_mass * 1000 / titer + + if abs(delta_caproic) > 1e-6 and total_other_base_yield > 1e-6: + adjustment_factor = max(0.1, (total_other_base_yield - delta_caproic) / total_other_base_yield) + for compound, base_yield in base_yields.items(): + if compound == 'CaproicAcid': + broth.imass[compound] = target_caproic_mass + else: + broth.imass[compound] = microalgae_mass * base_yield * adjustment_factor + else: + for compound, base_yield in base_yields.items(): + if compound == 'CaproicAcid': + broth.imass[compound] = target_caproic_mass + else: + broth.imass[compound] = microalgae_mass * base_yield + + current_volume = broth.F_vol + if current_volume > 0 and titer >= 5.0: + water_adjustment = max(0, target_volume_L - current_volume) + max_water = microalgae_mass * 10 + water_adjustment = min(water_adjustment, max_water) + broth.imass['Water'] += water_adjustment + gas.copy_like(n2) + gas.imass['H2'] = microalgae_mass * 0.1 + + def _design(self): + Design = self.design_results + Design['Fermenter size'] = self.outs[0].F_mass * self.tau + Design['Broth flow rate'] = self.outs[0].F_mass + duty = 50 * self.F_mass_in # Heat duty + self.add_heat_utility(duty, self.T) + +# %% +# ========================= +# Solid-Liquid Separation +# ========================= +@cost(basis='Flow rate', ID='Filtrate tank', units='kg/hr', + cost=103000, S=31815, CE=CEPCI[2010], n=0.7, BM=2.0) +@cost(basis='Flow rate', ID='Flitrate tank agitator', units='kg/hr', + kW=5.59275, cost=26000, S=337439, CE=CEPCI[2009], n=0.5, BM=1.5) +# ref from Biorefineries.HP +class SolidLiquidSeparation(Unit): + _N_ins = 1 + _N_outs = 2 + def _run(self): + feed = self.ins[0] + filtrate, residue = self.outs + residue.empty() + insolubles = ('Protein', 'Lipid', 'Ash', 'TriOlein', 'AlphaAmylase', 'GlucoAmylase', 'Yeast') + for comp in feed.chemicals.IDs: + if comp in insolubles: + residue.imass[comp] = feed.imass[comp] + filtrate.imass[comp] = 0 + else: + filtrate.imass[comp] = feed.imass[comp] + residue.imass[comp] = 0 + + filtrate.T = feed.T + filtrate.P = feed.P + residue.T = feed.T + residue.P = feed.P + + +# %% +# ========================= +# Other Units +# ========================= +@cost(basis='Flow rate', ID='NeutralizationTank', units='kg/hr', + # Size basis on the total flow, not just ammonia, + # thus assuming difference caused by MWs of NH3 and NH4OH is negligible + cost=5000, S=157478, CE=CEPCI[2009], n=0.5, BM=1.5) +class NeutralizationTank(Unit): + _N_ins = 2 + _N_outs = 1 + def _run(self): + acid_stream, ammonium_hydroxide = self.ins + self.outs[0].empty() + self.outs[0].mix_from([acid_stream, ammonium_hydroxide]) + neutralization_rxn = rxn.Rxn('2 NH4OH + H2SO4 -> AmmoniumSulfate + 2 H2O', 'NH4OH', 1) + neutralization_rxn(self.outs[0]) + + + + +@cost('Flow rate', 'Anaerobic Digestion Reactor', units= 'kg/hr', kW=100, S=339151, + cost=8000000, CE=CEPCI[2010], n=0.7, BM=1.8) +class AnaerobicDigestion(bst.Unit): + + _N_ins = 1 + _N_outs = 3 + + degrade_IDs = [ + 'AceticAcid', 'PropionicAcid', 'ButyricAcid', 'ValericAcid', + 'CaproicAcid', 'HeptanoicAcid', 'CaprylicAcid', + 'Ethanol', 'Butanol', 'Glucose', 'Protein', + 'AlphaAmylase', 'GlucoAmylase', 'TriOlein', 'WWTsludge' + ] + inorganic_IDs = ['Ash', 'H2SO4', 'Na2SO4', 'AmmoniumSulfate', 'H2O'] + + def __init__(self, ID='', ins=None, outs=(), *, tau=15*24, microalgae_mass: float|None=None, **kwargs): + super().__init__(ID, ins, outs, **kwargs) + self.tau = tau + self.microalgae_mass = microalgae_mass + self.T = 37 + 273.15 + + V_wf = 0.8 + + N_reactors = 0 + + N_transfer_pumps = 1 + + N_recirculation_pumps = 0 + + def _run(self): + feed = self.ins[0] + biogas, waste, sludge = self.outs + if self.microalgae_mass is None: + raise ValueError('microalgae_mass must be supplied when instantiating AnaerobicDigestion.') + microalgae_mass = self.microalgae_mass + H2_mol = microalgae_mass * 0.015 + CH4_mol = microalgae_mass * 0.015 + biogas.empty() + biogas.imol['H2'] = H2_mol + biogas.imol['CH4'] = CH4_mol + biogas.phase = 'g' + biogas.T = self.T + waste.empty() + sludge.empty() + for ID in self.degrade_IDs: + amount = feed.imass[ID] + waste.imass[ID] = amount * 0.05 # 5% into wastewater + sludge.imass[ID] = amount * 0.25 # 25% into sludge + for ID in self.inorganic_IDs: + amount = feed.imass[ID] + waste.imass[ID] = amount * 0.95 + sludge.imass[ID] = amount * 0.05 + sludge.imass['H2O'] = feed.imass['H2O']*0.3 + waste.imass['H2O'] = feed.imass['H2O']*0.7 + sludge.imass['Microalgae'] = 0 + waste.phase = 'l' + waste.T = sludge.T = self.T + + def _design(self): + Design = self.design_results + Design['Flow rate'] = sum([s.F_mass for s in self.outs]) # kg/hr + duty = 800 * self.microalgae_mass # Heat duty + self.add_heat_utility(duty, self.T) + + + + + + diff --git a/biorefineries/microalgae/utils.py b/biorefineries/microalgae/utils.py new file mode 100644 index 00000000..6fb4e636 --- /dev/null +++ b/biorefineries/microalgae/utils.py @@ -0,0 +1,117 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Sat July 07 13:15:00 2025 + +Microalgae biorefinery to produce medium chain fatty acids +by anaerobic fermentation without external electron donor addition + +References +---------- +[1] BioSTEAM Documentation: + https://biosteam.readthedocs.io/en/latest/API/thermosteam/Chemicals.html +[2] Cortes-Peña et al., BioSTEAM: A Fast and Flexible Platform for the Design, + Simulation, and Techno-Economic Analysis of Biorefineries under Uncertainty. + ACS Sustainable Chem. Eng. 2020, 8 (8), 3302–3310. +[3] 3-Hydroxypropionic acid biorefineries project: + https://github.com/BioSTEAMDevelopmentGroup/Bioindustrial-Park/tree/master/biorefineries/HP + +@author: Xingdong Shi +@version: 0.0.1 +""" + +import numpy as np +import pandas as pd +import biosteam as bst +import thermosteam as tmo +from ._chemicals import chems + +_kg_per_ton = 907.18474 +_lb_per_kg = 2.20462 +# Baseline feedstock flow rate +def get_feedstock_flow(dry_composition, moisture_content, dry_flow): + dry_array = chems.kwarray(dry_composition) + wet_flow = dry_flow / (1-moisture_content) + moisture_array = chems.kwarray(dict(Water=moisture_content)) + feedstock_flow = wet_flow * (dry_array*(1-moisture_content)+moisture_array) + return feedstock_flow + +# !!! This is the dry composition of microalgae; update to the composition you need +dry_composition = dict( + Carbohydrate=0.234, Lipid=0.059, Protein=0.5369, Ash=0.0512, Extract=0.1189) + +moisture_content = 0.05 #!!! This is the moisture content of preprocessed microalgae; update to the moisture content you need +dry_feedstock_flow = 2205. * _kg_per_ton / 24. # !!! Changing this will alter the total flow without altering the composition + +baseline_feedflow = get_feedstock_flow(dry_composition, moisture_content, + dry_feedstock_flow) + +# Price dictionary for all relevant chemicals and products +price = { + 'GlucoAmylase': 6.16, # $/kg (NREL 2018) + 'AlphaAmylase': 6.16, # $/kg (NREL 2018) + 'Yeast': 2, # $/kg Alibaba + 'Ethanol': 1.78, # $/kg https://catcost.chemcatbio.org/materials-library + 'OleylAlcohol': 1.00, # $/lb https://web.archive.org/web/20161125084558/http://www.icis.com:80/chemicals/channel-info-chemicals-a-z/ + 'AceticAcid': 1.31, # $/kg https://catcost.chemcatbio.org/materials-library + 'ButyricAcid': 1.72, # $/kg https://www.imarcgroup.com/butyric-acid-pricing-report? + 'CaproicAcid': 2.89, # $/kg Increasing the economic value of lignocellulosic stillage through medium-chain fatty acid production + 'Butanol': 2.27, # $/kg https://catcost.chemcatbio.org/materials-library + #'Butanol': 0.7, # $/lb https://web.archive.org/web/20161125084558/http://www.icis.com:80/chemicals/channel-info-chemicals-a-z/ + 'HeptanoicAcid': 3.17, # $/kg https://www.expertmarketresearch.com/prefeasibility-reports/heptanoic-acid-manufacturing-plant-project-report + 'CaprylicAcid': 5.07, # $/kg Increasing the economic value of lignocellulosic stillage through medium-chain fatty acid production + 'SulfuricAcid': 0.11, # $/kg https://catcost.chemcatbio.org/materials-library + 'AmmoniumHydroxide': 0.204, # BDO, succinic, HP program + 'NaOH': 1.01, #https://catcost.chemcatbio.org/materials-library + # Baseline from Davis et al., 2018, lower bound is 2015-2019 average of + # hydrate lime in $/ton at plant from Mineral Commodity Summaries 2020. + # 2015: 146.40 * (1.114/1.100) / 907.18474 = 0.163 + # 2016: 145.50 / 907.18474 = 0.160 + # 2017: 147.10 * (1.114/1.134) / 907.18474 = 0.159 + # 2018: 151.50 * (1.114/1.157) / 907.18474 = 0.161 + # 2019: 151.00 * (1.114/1.185) / 907.18474 = 0.156 + # (0.163+0.160+0.159+0.161+0.156) / 5 = 0.160 + # Upper bound is +10% from baseline = 0.1189 * _lb_per_kg * 1.1 = 0.288 + 'Lime': 0.1189 * _lb_per_kg +} + +# Chemical Engineering Plant Cost Index from Chemical Engineering Magzine +# (https://www.chemengonline.com/the-magazine/) +CEPCI = {1997: 386.5, + 1998: 389.5, + 2007: 525.4, + 2009: 521.9, + 2010: 550.8, + 2011: 585.7, + 2012: 584.6, + 2013: 567.3, + 2014: 576.1, + 2016: 541.7, + 2017: 567.5, + 2018: 603.1, + 2019: 607.5, + 2020: 596.2, + 2021: 708.8, + 2022: 816.0, + 2023: 797.9} + +_lps = bst.HeatUtility.get_heating_agent('low_pressure_steam') +_mps = bst.HeatUtility.get_heating_agent('medium_pressure_steam') +_hps = bst.HeatUtility.get_heating_agent('high_pressure_steam') +_mps.T = 233 + 273.15 +_hps.T = 266 + 273.15 + +_cooling = bst.HeatUtility.get_cooling_agent('cooling_water') +_chilled = bst.HeatUtility.get_cooling_agent('chilled_water') +_cooling.regeneration_price = 0 +_cooling.T = 28 + 273.15 +_cooling.T_limit = _cooling.T + 9 + +# Side steam in CHP not a heat utility, thus will cause problem in TEA utility +# cost calculation if price not set to 0 here, costs for regeneration of heating +# and cooling utilities will be considered as CAPEX and OPEX of CHP and CT, respectively +for i in (_lps, _mps, _hps, _cooling, _chilled): + i.heat_transfer_price = i.regeneration_price = 0 + + + diff --git a/biorefineries/microalgae/vis.py b/biorefineries/microalgae/vis.py new file mode 100644 index 00000000..889e356b --- /dev/null +++ b/biorefineries/microalgae/vis.py @@ -0,0 +1,797 @@ +""" +Created on Thu Jun 20 15:30:00 2024 + +Visualization module for Microalgae biorefinery to produce medium chain fatty acids +This script creates technical, economic and environmental figures for the system + +@author: Xingdong Shi +""" + +import os +import numpy as np +import pandas as pd +import matplotlib.pyplot as plt + +import biosteam as bst +import thermosteam as tmo +from biosteam.evaluation import Model, Metric +from chaospy import distributions as shape + +from .lca import create_microalgae_lca +from .system import create_microalgae_MCCA_production_sys, compute_labor_cost +from .tea import create_tea + +# Set style for publication quality plots +plt.rcParams['font.family'] = 'sans-serif' +plt.rcParams['font.sans-serif'] = ['Arial'] +plt.rcParams['font.size'] = 8 +plt.rcParams['axes.linewidth'] = 0.5 +plt.rcParams['axes.titlesize'] = 10 +plt.rcParams['axes.labelsize'] = 9 +plt.rcParams['xtick.labelsize'] = 8 +plt.rcParams['ytick.labelsize'] = 8 +plt.rcParams['legend.fontsize'] = 7 + +def calculate_breakdown_data(tea_obj, system): + """从TEA对象和系统中提取实际的成本分布数据""" + + try: + # 总投资成本 + total_TCI = tea_obj.TCI + total_VOC = tea_obj.VOC + total_FOC = tea_obj.FOC + + # 按类别分组单元 + feedstock_units = [] + conversion_units = [] + separation_units = [] + wastewater_units = [] + facilities_units = [] + + # 分类系统中的单元 + for unit in system.units: + unit_name = unit.__class__.__name__ + unit_id = unit.ID + + if 'U1' in unit_id or 'feed' in unit_id.lower(): + feedstock_units.append(unit) + elif 'R3' in unit_id or 'ferment' in unit_id.lower(): + conversion_units.append(unit) + elif any(x in unit_id for x in ['S3', 'D3', 'F3', 'M3']): + separation_units.append(unit) + elif 'waste' in unit_id.lower() or 'W' in unit_id: + wastewater_units.append(unit) + else: + facilities_units.append(unit) + + # 计算各类别的设备成本 + feedstock_cost = sum(getattr(unit, 'purchase_cost', 0) for unit in feedstock_units) + conversion_cost = sum(getattr(unit, 'purchase_cost', 0) for unit in conversion_units) + separation_cost = sum(getattr(unit, 'purchase_cost', 0) for unit in separation_units) + wastewater_cost = sum(getattr(unit, 'purchase_cost', 0) for unit in wastewater_units) + facilities_cost = sum(getattr(unit, 'purchase_cost', 0) for unit in facilities_units) + + # 归一化为百分比 + total_equipment_cost = feedstock_cost + conversion_cost + separation_cost + wastewater_cost + facilities_cost + if total_equipment_cost > 0: + equipment_costs = { + 'Feedstock': feedstock_cost / total_equipment_cost * 100, + 'Conversion': conversion_cost / total_equipment_cost * 100, + 'Separation': separation_cost / total_equipment_cost * 100, + 'Wastewater': wastewater_cost / total_equipment_cost * 100, + 'Cooling utility facilities': 0, + 'Boiler & turbogenerator': 0, + 'Other facilities': facilities_cost / total_equipment_cost * 100, + 'Heat exchanger network': 0, + 'Natural gas (steam)': 0, + 'Natural gas (drying)': 0, + 'Fixed operating costs': 0, + 'Credits': 0 + } + else: + # 回退到估算值 + equipment_costs = {k: np.nan for k in [ + 'Feedstock','Conversion','Separation','Wastewater', + 'Cooling utility facilities','Boiler & turbogenerator','Other facilities', + 'Heat exchanger network','Natural gas (steam)','Natural gas (drying)', + 'Fixed operating costs','Credits']} + + # 计算公用设施消耗 + total_cooling = total_heating = total_electricity = 0 + + for unit in system.units: + # 冷却负荷 + if hasattr(unit, 'heat_utilities'): + for hu in unit.heat_utilities: + if hu.duty < 0: + total_cooling += abs(hu.duty) + else: + total_heating += hu.duty + + # 电力消耗 + if hasattr(unit, 'power_utility') and unit.power_utility: + total_electricity += abs(unit.power_utility.power) + + # 按单元类型分配公用设施负荷 (简化分配) + cooling_duty = {k: np.nan for k in equipment_costs} + + heating_duty = {k: np.nan for k in equipment_costs} + + electricity = {k: np.nan for k in equipment_costs} + + # 运营成本分配 (基于FOC和VOC) + operating_cost = {k: np.nan for k in equipment_costs} + + # 输出实际数值供参考 + print(f"实际TEA数据:") + print(f" 总投资成本 (TCI): ${total_TCI/1e6:.1f} M") + print(f" 年度变动成本 (VOC): ${total_VOC/1e6:.1f} M/y") + print(f" 年度固定成本 (FOC): ${total_FOC/1e6:.1f} M/y") + print(f" 总加热负荷: {total_heating/1e6:.1f} MW") + print(f" 总冷却负荷: {total_cooling/1e6:.1f} MW") + print(f" 总电力消耗: {total_electricity/1e3:.1f} MW") + + except Exception as e: + print(f"获取实际成本数据时出错: {str(e)}") + + return { + 'Installed\nequipment\ncost': equipment_costs, + 'Cooling\nduty': cooling_duty, + 'Heating\nduty': heating_duty, + 'Electricity\nconsumption': electricity, + 'Operating\ncost': operating_cost + } + + +def create_plots_from_monte_carlo(save_path=None): + """Create and save publication-quality figures for microalgae MCCA system using Monte Carlo results.""" + + + try: + bst.main_flowsheet.clear() + bst.main_flowsheet.set_flowsheet('_') + + from microalgae._chemicals import chems + tmo.settings.set_thermo(chems) + except: + pass + + # 检查results目录是否存在，如果不存在则创建 + os.makedirs('results', exist_ok=True) + + # 定义结果文件路径 + mc_file_path = 'results/monte_carlo_results.csv' + + # 移除所有缓存检查，直接运行新模拟 + print("运行新的Monte Carlo模拟...") + # 初始化系统 + microalgae_sys = create_microalgae_MCCA_production_sys() + microalgae_sys.simulate() + + # 创建TEA + u = microalgae_sys.flowsheet.unit + s = microalgae_sys.flowsheet.stream + dry_tpd = u.U101.ins[0].F_mass * 24 / 1000 # kg/h -> t/d + microalgae_tea = create_tea( + system=microalgae_sys, + IRR=0.10, + duration=(2024, 2045), + depreciation='MACRS7', + income_tax=0.21, + operating_days=330, + lang_factor=None, + construction_schedule=(0.08, 0.60, 0.32), + startup_months=3, + startup_FOCfrac=1, + startup_salesfrac=0.5, + startup_VOCfrac=0.75, + WC_over_FCI=0.05, + finance_interest=0.08, + finance_years=10, + finance_fraction=0.4, + OSBL_units=(u.CT, u.CWP, u.ADP, u.PWC), + warehouse=0.04, + site_development=0.09, + additional_piping=0.045, + proratable_costs=0.10, + field_expenses=0.10, + construction=0.20, + contingency=0.10, + other_indirect_costs=0.10, + labor_cost=max(0.5e6, compute_labor_cost(dry_tpd)), + labor_burden=0.90, + property_insurance=0.007, + maintenance=0.03, + boiler_turbogenerator=None, + steam_power_depreciation='MACRS20' + ) + + # 设置主产品和相关变量，用于后续LCA计算 + main_product = s.caproic_acid_product + main_product_chemical_IDs = ['CaproicAcid'] + + # 尝试找到锅炉单元 + boiler = None + for unit in microalgae_sys.units: + if hasattr(unit, 'heat_utilities') and unit.heat_utilities: + boiler = unit + break + + # 如果没有找到锅炉，使用热交换器 + if boiler is None: + for unit in microalgae_sys.units: + if unit.__class__.__name__ == 'HXutility': + boiler = unit + break + + # 创建评估模型 + model = Model(microalgae_sys) + + # 定义setter函数 + def set_microalgae_price(price): + u.U101.ins[0].price = price + + def set_caproic_acid_price(price): + s.caproic_acid_product.price = price + + def set_caproic_acid_yield_factor(factor): + u.R301.caproic_acid_yield_factor = factor + + def set_operating_days(days): + microalgae_tea.operating_days = days + + def set_maintenance_factor(factor): + microalgae_tea.maintenance = factor + + # 定义指标 + model.metrics = [ + Metric('ROI', lambda: microalgae_tea.ROI * 100, '%'), + Metric('MSP', lambda: microalgae_tea.solve_price(s.caproic_acid_product), '$/kg'), + Metric('TCI', lambda: microalgae_tea.TCI, '$'), + Metric('VOC', lambda: microalgae_tea.VOC, '$'), + Metric('FOC', lambda: microalgae_tea.FOC, '$'), + Metric('GWP', lambda: create_microalgae_lca(system=microalgae_sys, + main_product=main_product, + main_product_chemical_IDs=main_product_chemical_IDs, + boiler=boiler).GWP, 'kg CO2 eq./kg') + ] + + # 注册参数 + model.parameter(setter=set_microalgae_price, element=u.U101.ins[0], + name='Microalgae Price', kind='isolated', + units='$/kg', distribution=shape.Triangle(-0.2, -0.1, 0.1)) + + model.parameter(setter=set_caproic_acid_price, element=s.caproic_acid_product, + name='Caproic Acid Price', kind='isolated', + units='$/kg', distribution=shape.Triangle(2.0, 2.89, 4.0)) + + model.parameter(setter=set_caproic_acid_yield_factor, element=u.R301, + name='Caproic Acid Yield Factor', kind='isolated', + distribution=shape.Triangle(0.7, 1.0, 1.3)) + + model.parameter(setter=set_operating_days, element=microalgae_tea, + name='Operating Days', kind='isolated', + distribution=shape.Triangle(300, 330, 350)) + + model.parameter(setter=set_maintenance_factor, element=microalgae_tea, + name='Maintenance Factor', kind='isolated', + distribution=shape.Triangle(0.02, 0.03, 0.05)) + + # 运行Monte Carlo模拟 + n_samples = 10 + + # 获取参数分布 + param_distributions = {} + for param in model.parameters: + if hasattr(param, 'distribution'): + param_distributions[param.index] = param.distribution + + # 获取指标名称映射 + metric_names = {} + for i, metric in enumerate(model.metrics): + metric_names[i] = metric.name + + # 获取基准样本作为模板 + baseline = model.get_baseline_sample() + results = [] + + # 生成样本并评估 + for i in range(n_samples): + if i % 100 == 0 and i > 0: + print(f"Completed {i}/{n_samples} samples") + + # 复制基准样本 + sample = baseline.copy() + + # 为每个参数生成随机值 + for param_name, distribution in param_distributions.items(): + # 简化的拉丁超立方抽样 + u_val = (i + 0.5) / n_samples + if hasattr(distribution, 'ppf'): + sample[param_name] = distribution.ppf(u_val) + else: + sample[param_name] = distribution.sample() + + # 评估样本 + try: + result_series = model(sample) + + # 创建包含参数和指标值的字典 + result_dict = {} + + # 添加参数值 + for param_name in param_distributions.keys(): + result_dict[param_name] = sample[param_name] + + # 添加指标值 + for i, name in metric_names.items(): + result_dict[name] = result_series.iloc[i] + + # 保障 GWP 指标存在（某些情况下 Metric 计算可能失败） + if 'GWP' not in result_dict or np.isnan(result_dict['GWP']): + try: + lca_tmp = create_microalgae_lca(system=microalgae_sys, + main_product=main_product, + main_product_chemical_IDs=main_product_chemical_IDs, + boiler=boiler) + result_dict['GWP'] = lca_tmp.GWP + except Exception as _: + result_dict['GWP'] = np.nan + + results.append(result_dict) + except Exception as e: + print(f"Sample {i} evaluation failed: {str(e)}") + + # 组合结果 + monte_carlo_results = pd.DataFrame(results) + + # 保存结果 + monte_carlo_results.to_csv(mc_file_path, index=False) + + # ========================================================================= + # FIGURE A: MFSP comparison across different scales (Monte Carlo based) + # ========================================================================= + fig = plt.figure(figsize=(7, 8)) + + # A: MFSP across different scales + ax1 = plt.subplot2grid((3, 2), (0, 0)) + + # 使用 Monte Carlo 样本结果绘制 MFSP 分布 (均值 ± 标准差) + if 'MSP' in monte_carlo_results.columns: + mfsp_series = monte_carlo_results['MSP'].dropna() + mean_mfsp = mfsp_series.mean() + std_mfsp = mfsp_series.std() + ax1.bar(['Monte Carlo'], [mean_mfsp], yerr=[std_mfsp], color='#A52A2A', alpha=0.7, width=0.5) + # 价格范围背景 + ax1.axhspan(2.5, 3.5, alpha=0.3, color='gray', label='Bio-based MFSP range') + ax1.axhspan(1.0, 2.5, alpha=0.2, color='gray', label='Market price range') + ax1.set_ylim(0, max(5.0, mean_mfsp*1.5)) + else: + ax1.text(0.5, 0.5, 'No MSP data', ha='center', va='center') + + ax1.set_ylabel('MFSP [$kg^{-1}$]') + ax1.set_title('MFSP (Monte Carlo mean ± std)') + ax1.set_xticks([]) + + legend_elements = [ + plt.Rectangle((0,0), 1, 1, fc='gray', alpha=0.3, label='Bio-based MFSP range'), + plt.Rectangle((0,0), 1, 1, fc='gray', alpha=0.2, label='Market price range') + ] + ax1.legend(handles=legend_elements, loc='upper left', bbox_to_anchor=(1.05, 1), fontsize=7) + + # Cost and utility breakdown chart - 使用计算数据 + ax2 = plt.subplot2grid((3, 2), (0, 1)) + + # 通过计算获取实际的成本和公用设施数据 + # 获取分类数据 + breakdown_data = calculate_breakdown_data(microalgae_tea, microalgae_sys) + + # 定义类别 + categories = ['Installed\nequipment\ncost', 'Cooling\nduty', 'Heating\nduty', 'Electricity\nconsumption', 'Operating\ncost'] + + # 创建数据字典 + data = {} + for key in list(breakdown_data.values())[0].keys(): # 使用第一个字典的键 + data[key] = [] + for category in categories: + data[key].append(breakdown_data[category][key]) + + # 创建DataFrame用于绘图 + df = pd.DataFrame(data, index=categories) + + # 保存成本与公用设施百分比分布到 CSV + try: + df.to_csv(os.path.join('results', 'cost_utility_breakdown_percent.csv')) + except Exception as _: + pass + + # 绘制堆叠条形图 + bottom = np.zeros(len(categories)) + colors = ['#FADBD8', '#A52A2A', '#641E16', '#FFFFFF', '#EAECEE', '#D5D8DC', + '#85929E', '#2C3E50', '#1C2833', '#17202A', '#7D3C98', '#FFFFFF'] + hatches = ['', '', '', '', '', '', '', '', '', '', '', '////'] + + for i, (col, color, hatch) in enumerate(zip(df.columns, colors, hatches)): + ax2.bar(categories, df[col], bottom=bottom, label=col, color=color, hatch=hatch) + bottom += df[col] + + # 计算每个类别的总和和单位 - 使用实际数据 + try: + # 从TEA获取实际数值 + actual_TCI = microalgae_tea.TCI / 1e6 # M$ + actual_heating = sum(abs(hu.duty) for unit in microalgae_sys.units + if hasattr(unit, 'heat_utilities') + for hu in unit.heat_utilities if hu.duty > 0) / 1e6 # MW + actual_cooling = sum(abs(hu.duty) for unit in microalgae_sys.units + if hasattr(unit, 'heat_utilities') + for hu in unit.heat_utilities if hu.duty < 0) / 1e6 # MW + actual_electricity = sum(abs(unit.power_utility.power) for unit in microalgae_sys.units + if hasattr(unit, 'power_utility') and unit.power_utility) / 1e3 # MW + actual_operating = (microalgae_tea.VOC + microalgae_tea.FOC) / 1e6 # M$/y + + sum_values = [actual_TCI, actual_cooling, actual_heating, actual_electricity, actual_operating] + except Exception as e: + print(f"[WARNING] Failed to fetch aggregate totals: {e}") + sum_values = [np.nan]*5 + + units = ['[$10^6$ $]', '[MW]', '[MW]', '[MW]', '[$10^6$ $/y]'] + + # 添加总和值和单位 + for i, category in enumerate(categories): + if i < len(sum_values) and not np.isnan(sum_values[i]): + plt.text(i, 105, f'sum: {sum_values[i]:.1f}', ha='center', fontsize=6) + plt.text(i, 112, units[i], ha='center', fontsize=6) + + ax2.set_ylim(-40, 120) + ax2.set_ylabel('Cost and Utility Breakdown [%]') + # 将图例移到图右侧 + ax2.legend(loc='upper left', bbox_to_anchor=(1.05, 1.0), fontsize=5) + + # ========================================================================= + # FIGURE B: GWP comparison and breakdown (使用LCA计算) + # ========================================================================= + + # B: GWP across different scales - 使用LCA计算的实际数据 + ax3 = plt.subplot2grid((3, 2), (1, 0)) + + # 使用与MFSP相同的条件计算GWP值 + # 这里复用之前设置的conditions + if 'GWP' in monte_carlo_results.columns: + gwp_series = monte_carlo_results['GWP'].dropna() + mean_gwp = gwp_series.mean() + std_gwp = gwp_series.std() + ax3.bar(['Monte Carlo'], [mean_gwp], yerr=[std_gwp], color='#A52A2A', alpha=0.7, width=0.5) + ax3.set_ylim(0, max(5.0, mean_gwp*1.3)) + else: + ax3.text(0.5,0.5,'No GWP data',ha='center',va='center') + ax3.set_ylim(0,5) + + # 添加参考线 - 行业基准值 + f1_value = 3.1 # 行业平均值示例 + f2_value = 1.7 # 行业最佳实践示例 + ax3.axhline(y=f1_value, linestyle='--', color='black', linewidth=0.5) # f1 + ax3.axhline(y=f2_value, linestyle='--', color='black', linewidth=0.5) # f2 + + # 文字标注 + ax3.text(2.1, f1_value, '$f_1$', fontsize=8) + ax3.text(2.1, f2_value, '$f_2$', fontsize=8) + + ax3.set_ylabel('GWP$_{100}$ [kg CO$_2$ eq. kg$^{-1}$]') + + # GWP breakdown chart - 使用LCA的实际数据 + ax4 = plt.subplot2grid((3, 2), (1, 1)) + + # 从LCA获取GWP明细 + # 重新计算当前状态的LCA + lca = create_microalgae_lca( + system=microalgae_sys, + main_product=main_product, + main_product_chemical_IDs=main_product_chemical_IDs, + boiler=boiler + ) + + # 将数据转换为百分比的字典 + gwp_total = abs(lca.GWP) if lca.GWP != 0 else 1.0 # 避免除零错误 + gwp_data = { + 'Feedstock': [max(0, lca.feedstock_GWP/gwp_total*100)], + 'Material Inputs': [max(0, lca.material_GWP/gwp_total*100)], + 'Electricity': [max(0, lca.net_electricity_GWP/gwp_total*100)], + 'Direct Emissions': [max(0, lca.direct_non_biogenic_emissions_GWP/gwp_total*100)], + } + + # 如果有负值(碳捕获)，添加Credits类别 + if lca.feedstock_GWP < 0: + gwp_data['Credits'] = [min(0, lca.feedstock_GWP/gwp_total*100)] # 负值作为贷项 + else: + gwp_data['Credits'] = [0] + + # 创建DataFrame以进行绘图 + gwp_df = pd.DataFrame(gwp_data) + + # 保存 GWP breakdown 百分比到 CSV + try: + gwp_df.to_csv(os.path.join('results', 'gwp_breakdown_percent.csv'), index=False) + except Exception as _: + pass + + # 绘制堆叠条形图 - GWP + gwp_bottom = np.zeros(1) + gwp_colors = ['#FADBD8', '#A52A2A', '#1C2833', '#85929E', '#FFFFFF'] + gwp_hatches = ['', '', '', '/////', '////'] + + for i, (col, color, hatch) in enumerate(zip(gwp_df.columns, gwp_colors, gwp_hatches)): + ax4.bar(['GWP$_{100}$ Breakdown'], gwp_df[col], bottom=gwp_bottom, label=col, color=color, hatch=hatch) + gwp_bottom += gwp_df[col] + + ax4.set_ylim(-40, 100) + ax4.set_ylabel('GWP$_{100}$ Breakdown [%]') + # 将图例移到图右侧 + ax4.legend(loc='upper left', bbox_to_anchor=(1.05, 1.0), fontsize=5) + + # 添加脚注文本 - 基于实际LCA数据 + footnote = "*Feedstock (microalgae) farming, harvesting,\n"\ + "transportation, storage, and handling. Credit for\n"\ + f"fixed carbon: {lca.feedstock_GWP:.2f} kg CO2 eq./kg.\n"\ + "Sum of direct biogenic emissions: "\ + f"{lca.biogenic_emissions_GWP:.2f} kg CO2 eq./kg." + plt.figtext(0.1, 0.02, footnote, fontsize=7, ha='left') + + # ========================================================================= + # FIGURES C-F: Contour plots for yield vs titer (使用实际计算数据) + # ========================================================================= + + # 定义产率和浓度的网格点 - 与组合图保持一致 + yields = np.linspace(30, 80, 10) # % theoretical - 与组合图一致 + titers = np.linspace(30, 100, 10) # g/L - 与组合图一致 + + # 创建数据存储数组 + Z_neutral_mfsp = np.zeros((len(titers), len(yields))) + Z_lowph_mfsp = np.zeros((len(titers), len(yields))) + Z_neutral_gwp = np.zeros((len(titers), len(yields))) + Z_lowph_gwp = np.zeros((len(titers), len(yields))) + + # 创建网格 - 与组合图一致 + Y, T = np.meshgrid(yields, titers) + + # 保存原始参数 + original_yield_factor = u.R301.caproic_acid_yield_factor if hasattr(u.R301, 'caproic_acid_yield_factor') else 1.0 + original_titer = u.R301.titer if hasattr(u.R301, 'titer') else 60.0 + + print("开始计算等高线图数据...") + # 计算每个网格点的MFSP和GWP值 + for i, titer in enumerate(titers): + for j, yield_pct in enumerate(yields): + try: + # 改进产率因子转换公式 - 更合理地映射产率到产率因子 + yield_factor = 0.7 + (yield_pct - 30) / (80 - 30) * 0.6 # 映射到0.7-1.3范围 + + # 设置参数 + if hasattr(u.R301, 'caproic_acid_yield_factor'): + u.R301.caproic_acid_yield_factor = yield_factor + + # 设置滴度参数 + if hasattr(u.R301, 'titer'): + u.R301.titer = titer + + # 中性发酵条件 + microalgae_sys.simulate() + + # 创建LCA计算GWP + try: + neutral_lca = create_microalgae_lca( + system=microalgae_sys, + main_product=main_product, + main_product_chemical_IDs=main_product_chemical_IDs, + boiler=boiler + ) + neutral_gwp = neutral_lca.GWP + except Exception as e: + # 如果LCA计算失败，使用简单估算模型 + print(f"LCA计算失败：{str(e)}") + raise + + # 计算MFSP + mfsp = microalgae_tea.solve_price(s.caproic_acid_product) + Z_neutral_mfsp[i, j] = mfsp + Z_neutral_gwp[i, j] = neutral_gwp + + # 低pH发酵条件 (假设pH对产率有5%的影响) + if hasattr(u.R301, 'caproic_acid_yield_factor'): + u.R301.caproic_acid_yield_factor = yield_factor * 0.95 + + microalgae_sys.simulate() + + # 创建LCA计算GWP + try: + lowph_lca = create_microalgae_lca( + system=microalgae_sys, + main_product=main_product, + main_product_chemical_IDs=main_product_chemical_IDs, + boiler=boiler + ) + lowph_gwp = lowph_lca.GWP + except Exception as e: + # 如果LCA计算失败，使用简单估算模型 + print(f"低pH LCA计算失败：{str(e)}") + raise + + mfsp_lowph = microalgae_tea.solve_price(s.caproic_acid_product) + Z_lowph_mfsp[i, j] = mfsp_lowph + Z_lowph_gwp[i, j] = lowph_gwp + + # 仅显示部分计算结果以减少输出 + if (i == 0 or i == len(titers)-1) and (j == 0 or j == len(yields)-1): + print(f"计算完成: 浓度={titer:.1f}, 产率={yield_pct:.1f}%, 产率因子={yield_factor:.3f}, MFSP={mfsp:.2f}, GWP={neutral_gwp:.2f}") + + except Exception as e: + print(f"[WARNING] Convergence failure (T={titer:.1f}, Y={yield_pct:.1f}%): {e}") + Z_neutral_mfsp[i, j] = np.nan + Z_lowph_mfsp[i, j] = np.nan + Z_neutral_gwp[i, j] = np.nan + Z_lowph_gwp[i, j] = np.nan + + # 恢复原始参数 + if hasattr(u.R301, 'caproic_acid_yield_factor'): + u.R301.caproic_acid_yield_factor = original_yield_factor + if hasattr(u.R301, 'titer'): + u.R301.titer = original_titer + microalgae_sys.simulate() + + # 添加平滑处理并在此处保存等高线数据矩阵 + try: + from scipy.ndimage import gaussian_filter + Z_neutral_mfsp = gaussian_filter(Z_neutral_mfsp, sigma=0.7) + Z_lowph_mfsp = gaussian_filter(Z_lowph_mfsp, sigma=0.7) + Z_neutral_gwp = gaussian_filter(Z_neutral_gwp, sigma=0.7) + Z_lowph_gwp = gaussian_filter(Z_lowph_gwp, sigma=0.7) + print("数据平滑处理已应用") + except ImportError: + print("无法导入scipy.ndimage，跳过数据平滑处理") + + # 保存等高线数据到 NPZ + try: + np.savez(os.path.join('results', 'contour_data.npz'), + yields=yields, titers=titers, + Z_neutral_mfsp=Z_neutral_mfsp, Z_lowph_mfsp=Z_lowph_mfsp, + Z_neutral_gwp=Z_neutral_gwp, Z_lowph_gwp=Z_lowph_gwp) + except Exception as _: + pass + + # Plot C: Neutral fermentation MFSP + ax5 = plt.subplot2grid((3, 2), (2, 0)) + plt.title('Neutral Fermentation') + cs = ax5.contourf(Y, T, Z_neutral_mfsp, 30, cmap='viridis_r') + + # 计算合适的轮廓线级别 + try: + # 计算均匀分布的8个轮廓线级别，确保递增 + min_val = np.nanmin(Z_neutral_mfsp) + max_val = np.nanmax(Z_neutral_mfsp) + print(f"MFSP数据范围: {min_val:.3f} - {max_val:.3f}") + + if min_val >= max_val or np.isnan(min_val) or np.isnan(max_val): + # 如果数据范围有问题，使用默认值 + levels = [1.0, 1.5, 2.0, 2.5, 3.0] + print("使用默认级别") + else: + levels = np.linspace(min_val, max_val, 8) + # 确保级别是递增的且唯一 + levels = np.unique(np.sort(levels)) + if len(levels) < 3: # 如果唯一级别太少，使用默认值 + levels = np.linspace(min_val, max_val, 5) + levels = np.round(levels, 3) + + print(f"MFSP等高线级别: {levels}") + # 验证级别是递增的 + if len(levels) < 2 or not np.all(np.diff(levels) > 0): + levels = [1.0, 1.5, 2.0, 2.5, 3.0] + print("回退到默认级别") + except Exception as e: + print(f"级别计算错误: {e}") + levels = [1.0, 1.5, 2.0, 2.5, 3.0] + + contours = ax5.contour(Y, T, Z_neutral_mfsp, levels=levels, colors='black', linewidths=0.5) + ax5.clabel(contours, inline=True, fontsize=6, fmt='%.2f') + ax5.set_xlabel('Yield [%theoretical]') + ax5.set_ylabel('Titer [g L$^{-1}$]') + + # Add extra space on the right and place colorbar there + fig.subplots_adjust(right=0.85) + cbar = fig.colorbar(cs, ax=ax5, location='right', pad=0.02) + + # 标记基准案例点 + base_yield_pct = 50.0 # 基准产率百分比 + base_titer = original_titer # 基准浓度 + ax5.plot(base_yield_pct, base_titer, 'D', color='white', markersize=5, markeredgecolor='black', markeredgewidth=0.5) + + # Adjust layout + plt.tight_layout() + fig.subplots_adjust(hspace=0.4, wspace=0.3, left=0.1, right=0.9, top=0.95, bottom=0.1) + + # Save the figure if a path is provided + if save_path: + plt.savefig(save_path, dpi=300, bbox_inches='tight') + + # Create additional GWP contour plots + fig2 = plt.figure(figsize=(7, 3)) + + # Plot E: Neutral fermentation GWP + ax7 = plt.subplot(1, 2, 1) + plt.title('Neutral Fermentation') + cs = ax7.contourf(Y, T, Z_neutral_gwp, 30, cmap='viridis_r') + + # 计算合适的轮廓线级别 + try: + # 计算均匀分布的8个轮廓线级别，确保递增 + min_val = np.nanmin(Z_neutral_gwp) + max_val = np.nanmax(Z_neutral_gwp) + print(f"GWP数据范围: {min_val:.3f} - {max_val:.3f}") + + if min_val >= max_val or np.isnan(min_val) or np.isnan(max_val): + # 如果数据范围有问题，使用默认值 + levels = [1.5, 2.0, 2.5, 3.0, 3.5] + print("使用默认级别") + else: + levels = np.linspace(min_val, max_val, 8) + # 确保级别是递增的且唯一 + levels = np.unique(np.sort(levels)) + if len(levels) < 3: # 如果唯一级别太少，使用默认值 + levels = np.linspace(min_val, max_val, 5) + levels = np.round(levels, 3) + + print(f"GWP等高线级别: {levels}") + # 验证级别是递增的 + if len(levels) < 2 or not np.all(np.diff(levels) > 0): + levels = [1.5, 2.0, 2.5, 3.0, 3.5] + print("回退到默认级别") + except Exception as e: + print(f"级别计算错误: {e}") + levels = [1.5, 2.0, 2.5, 3.0, 3.5] + + contours = ax7.contour(Y, T, Z_neutral_gwp, levels=levels, colors='black', linewidths=0.5) + ax7.clabel(contours, inline=True, fontsize=6, fmt='%.2f') + ax7.set_xlabel('Yield [%theoretical]') + ax7.set_ylabel('Titer [g L$^{-1}$]') + + # 添加基准点 + ax7.plot(base_yield_pct, base_titer, 'D', color='white', markersize=5, markeredgecolor='black', markeredgewidth=0.5) + + + # Add colorbar outside of the two subplots + fig2.subplots_adjust(right=0.85) + cbar = fig2.colorbar(cs, ax=ax7, location='right', pad=0.02) + + # Adjust layout + fig2.subplots_adjust(wspace=0.3, left=0.1, right=0.9) + + # Save the second figure if a path is provided + if save_path: + plt.savefig(save_path.replace('.png', '_gwp.png'), dpi=300, bbox_inches='tight') + + # ------------------------------------------------------------------ + # Additionally split every Axes object in fig and fig2 into its own + # image file so that users can access each sub-plot independently. + # ------------------------------------------------------------------ + base_dir = os.path.dirname(save_path) or '.' + split_dir = os.path.join(base_dir, 'combined_subplots') + os.makedirs(split_dir, exist_ok=True) + + # Helper to save each axis + def _save_axes_individually(fig_obj, tag): + fig_obj.canvas.draw() # Ensure renderer exists + renderer = fig_obj.canvas.get_renderer() + for i, ax in enumerate(fig_obj.axes, start=1): + # Use axis title as part of filename when available, else index + title = ax.get_title() or f'{tag}_{i}' + safe_title = ''.join(c if c.isalnum() else '_' for c in title)[:40] + fname = f'{tag}_{safe_title or i}.png' + path = os.path.join(split_dir, fname) + bbox = ax.get_tightbbox(renderer).transformed(fig_obj.dpi_scale_trans.inverted()) + fig_obj.savefig(path, dpi=300, bbox_inches=bbox) + + _save_axes_individually(fig, 'fig1') + _save_axes_individually(fig2, 'fig2') + + return fig, fig2 + +create_plots_from_monte_carlo('microalgae_mcca_summary.png') \ No newline at end of file